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1 .  INTRODUCTION 

Given  a  positive  number  q«1  let  denote  the  (upper)  qth  quantile 
of  a  random  variable  X,  defined  by  P(X>Xq)=q.  In  a  previous  report  (Breiman, 
Stone  and  Gins  [2])  a  number  of  point  estimators  of  xq  based  on  a  random 
sample  of  size  n  from  the  distribution  of  X  were  studied.  The  main  purpose 
of  the  present  research  is  to  develop  and  study  the  performance  of  several 
confidence  interval  procedures  for  x^  based  on  the  sample,  especially  when 
q=l/n.  This  work  is  reported  on  in  Section  2.  In  Subsection  2.1.1,  modifi¬ 
cations  of  one  of  these  procedures  are  developed  for  hand! ing  grouped  data, 
which  Contain  many  ties,  satisfactorily.  A  detailed  discussion  of  the 
quadratic  tail  procedure,  described  in  Subsection  2.3,  is  given  in  Section  3. 
In  Section  4  the  behavior  of  the  exponential  tail  estimator  is  studied  when 
the  sample  data  is  stationary  but  dependent.  The  selection  of  procedures 
for  consideration  has  been  guided  by  the  dictum  of  DuMouchel  and  Olshen  [3] 
that  one  should  "let  the  tails  of  the  data  speak  for  themselves." 
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2.  CONFIDENCE  INTERVALS  FOR  EXTREME  QUANTILES 

Let  X  denote  a  random  variable  whose  (unknown)  underlying  distribution 
function  F  is  continuous  on  (-“,»)  and  strictly  increasing  on  an  interval 
I  which  contains  the  support  of  F  [i.e.,  which  is  such  that  P(XeI)=l].  Let 
S(x),  x>0,  be  the  tail  (survival)  probability  function  defined  by 
S(x)=P(X>x)=1-F(x) .  For  0<q<l  let  xq  denote  the  (upper)qth  quantile  of  X 
defined  by 


s(Xq)  =  1-F(xq)  =  P(X>xq)  =  q 

Let  X-j , . . .  »Xn  denote  a  random  sample  of  size  n  from  F  and  let 

X^i  j,... ,X^nj  denote  the  corresponding  order  statistics  defined  so  that 

X^-j  )>*  •  *>x(n)-  A  confidence  interval  procedure  for  xq  is  an  interval 

I  =  [Xq,xq],  where  Xq<.xq  with  both  x^  and  7q  being  functions  of  X^  j .  ,X^nj. 

The  length  of  |I|  of  the  interval  I  is  given  by  I  =  7  -x^. 

Let  Pp  and  Ep  denote  probabilities  and  expectations  when  F  is  the 

distribution  function  of  X.  Relevant  characteristics  of  the  confidence 

interval  procedure  include  Pp(xq<T),  Pp(xq>T),  Pp(xq£T)  =  Pp(xq<T)  +  Pp(x  >T), 
and  Ep | I | .  (Here  x  <T  and  x  >1  mean  respectively  that  xa<>L  and  xa>xQ.)  The 

goals  of  making  the  indicated  probabilities  and  expected  length  both  small 

over  a  wide  range  of  F’s  are  obviously  in  conflict  with  each  other. 

One  approach  to  obtaining  a  confidence  interval  procedure  is  to  start 

out  with  a  realistic  model  F(*;0),  9e®,  for  F,  express  xq  as  a  function  xq(9) 

of  the  unknown  parameter  8,  and  then  employ  a  classical  parametric  confidence 

interval  procedure  for  xq(9)-  Given  0<a<l,  suppose  T  is  (at  least  approximately) 

a  classical  100(l-a}%  confidence  interval  based  on  the  assumed  model;  i.e.,  that 
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Pe(yi)  i  a,  Be®  , 

where  PQ  -  Pp^..Qj.  Then  roughly  speaking,  Pp(xqtfT)  •  a  for  any  F  which 
can  be  globally  well  approximated  by  F(*;B)  for  some  6e®.  Otherwise 
Pp(xq^T)  can  be  substantially  larger  than  a.  This  is  true  in  particular 
for  F's  whose  central  portion  and  extreme  upper  tail  are  best  approximated 
by  distribution  functions  F( * ;0 )  with  significantly  different  values  of  9. 
Relative  to  such  F's  the  confidence  interval  procedure  is  not  robust  vis 
a  vis 


Pp(xq^I)  •  a  .  (2.1) 

There  are  two  ways  of  making  such  a  confidence  procedure  more  robust: 
(1)  Consider  a  higher  dimensional  model  F(*-,0,x),  (8,x)e(H)xT  (e.g.,  a  three- 
parameter  instead  of  a  two-parameter  lognormal  model).  A  wider  range  of 
distribution  functions  F  can  be  globally  well  approximated  by  a  distribution 
function  in  this  larger  model.  (2)  Base  the  confidence  interval  procedure 
on  only  the  upper  m  order  statistics  X^,...,X^  for  some  nKn,  rather  than 
on  all  the  original  data.  Both  (1)  and  (2)  lead  to  procedures  which  are  more 
robust  vis  a  vis  Eq.  (2.1).  Unfortunately  they  also  both  lead  to  significant 
increases  in  Ep  T  . 

Robustness  vis  a  vis  Eq.  (2.1)  here  is  with  respect  to  departures  of  F 
from  the  assumed  model.  Another  type  of  robustness,  when  F  belongs  to  the 
assumed  model,  is  with  respect  to  errors  in  measuring  or  recording  the  sample 
data;  this  type  of  robustness  will  not  be  considered  in  the  present  report. 
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[Note  that  (1)  and  (2)  above  can  lead  to  procedures  which  are  less  (not  more! ) 
robust  with  respect  to  measurement  and  recording  errors.] 

Understandability  and  ease  of  implementation  are  important  additional 
considerations  in  determining  which  procedures  to  use. 

So  far,  confidence  intervals  I=I(q)  have  been  considered  for  quantiles 
Xq,  0<q<l .  It  is  also  desirable  to  consider  confidence  intervals  7=J(x)  for 
tail  probabilities  S(x)=l-F(x),  -®<x<®.  There  is  a  natural  one-to-one 
correspondence  between  confidence  intervals  for  xq  and  these  for  l-F(x) 
given  by  7(x)={q:xeT(q)}  and  T(q)={x:qe7(x)}.  This  correspondence  preserves 
coverage  probabilities.  That  is, 

PF(YT(q))  =  PF(qeJ(xq)) 
and 

PF(l-F(x)eT(x))  =  PF(x1.F(x)eT0-F(x)  ))  if  0<F(x)<l  . 

Because  of  this  close  correspondence  between  confidence  intervals  for  quantiles 
and  those  for  tail  probabilities  it  was  decided  to  devote  the  present  research 
effort  exclusively  to  confidence  intervals  for  quantiles. 

Three  confidence  interval  procedures  will  be  described  in  Subsections  2.1 
through  2.3.  (The  procedure  described  in  Subsection  2.3  will  be  elaborated 
on  in  Section  3.)  A  Monte  Carlo  experiment  designed  to  compare  the  performance 
characteri sties  of  these  procedures  will  be  discussed  in  Subsection  2.4. 
Tentative  conclusions  drawn  from  the  results  of  this  experiment  and  suggestions 
for  further  work  are  presented  in  Subsection  2.5. 
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2.1  TWO-PARAMETER  EXPONENTIAL  PROCEDURE 


Consider  the  two-parameter  exponential  model 
F(x;x,a)  =  l-e~^x~T^a,  x>x, 

=  0,  x< T, 

where  Te(-«,o°)  is  a  location  parameter  and  a>0  is  a  scale  parameter. 
Correspondingly 

S(x;T,a)  =  e~^x~T^a,  x>r, 


and 


»  1, 


X<T  , 


xq  =  t+a  log(l/q),  0<q<l 

The  two-parameter  exponential  model  is  appealing  for  several  reasons. 

First,  it  is  very  simple  and  leads  to  the  above  simple  formula  for  xq. 
Second,  the  upper  tail  of  distributions  of  this  type  can  be  used  to  provide 
reasonably  accurate  approximations  to  the  upper  tail  of  a  number  of  commonly 
assumed  alternative  model s--Wei bull ,  gamma,  and  lognormal.  Third,  the  upper 
tail  of  distributions  of  this  type  is  realistic  in  many  applications  (see, 
e.g.,  Breiman,  Gins,  and  Stone  [1]). 
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Given  2<m<n,  let  denote  the  tipper  m  order  statistics 

based  on  a  random  sample  of  size  n  from  F(*;x,a).  The  joint  density  of  these 
random  variables  is  given  by 


fy  y  ( X-, , . . .  ,x  ;t »p) 

X(i),...,X(m)  1  m 


hi 

l(xrxm)/a  -m(xm-x)/a 


n:  .-m.  l~i  "m"  “  ’"'"m  l/'“  f,  "^VT^/a 

IrwriTT  a  e  1  e  |_ 


n-m 


for  x<x  <. . .<x, ,  while  fv  v  (x, , . . . ,x  :x,a)=0  otherwise.  The 

m  1  X(1 )’*■*’  X{m)  1  m 

maximum  likelihood  estimators  of  x  and  a  based  on  X,,»,...,X/  a  are  given  by 

(1)  (m; 


a=i^  CX(i)“X(m)] 


m-1 


and 


r  =  X(m) 


The  corresponding  maximum  likelihood  estimator  of  xq  is  given  by 


xq  =  T  ♦  a  log  I  =  X(m)  +  a  log  i 


This  estimator,  with  a  replaced  by 


m-1 


=  iry  £  Cx(i)-x(m)] 
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was  studied  in  detail  in  Breiman,  Stone,  and  Gins  [2],  where  it  was  referred 
to  as  the  exponential  tail  estimator.  Considering  its  simplicity,  it  is 
surprisingly  robust  to  moderate  departures  of  the  upper  tail  of  F  from  the 
assumed  exponential  form. 

A  confidence  interval  procedure  for  x^  will  now  be  obtained.  To  this 
end  set  »,-X(1)-*(b)  for  1  <i<m-l .  Then 


fy  y  x  (y-j  >  •  •  •  ) 

Y1 . Vl,x(m)  1  m 


n:  -m  -I,  va  -m<vT)/ar.  -<VT)/a] 

rara  ole  l-e 


for  xm>T  and  0<ynl_1<-..<y1.  while  ,Ym_,  ,*(nl)(yl . VrVT'a>'° 

otherwise.  Let  Zp...,Z  ^  be  independent  random  variables  having  the 

common  exponential  density  f,  defined  by 

Z1 


f7  (z)  =  \  e"z/a  ,  z>0 

Z1  a 

=  0  ,  z<0  , 

and  suppose  that  ,. . .  ,Zm_-|  ,X^  are  independent  random  variables.  Let 
j ,. . .  j  be  the  order  statistics  from  Z-j ,. . .  ,Zm_^ ,  defined  so  that 

Z^  1  j>. . . >Z(m_-|  ^ •  Then  Z^  j ,. . .  ,Z^m_i  j  ,X^mj  has  the  same  joint  distribution 
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as  Y-j »X(m)  •  Consequently  Y-|+- •  .+Ym_i  is  independent  of  and 
has  the  same  distribution  as 

Z(1 )+* " '+Z(m-1 )  =  Z]+- • ,+zm-i 

namely  the  gamma  distribution  having  density  g(*;m-l,l/a)  defined  by 


.m-2  -t/a 

g(t;m-l  ,1  /a)  =  -  ,  t>0 

am“  (m-1 )1 

=  0  ,  t<0 

Therefore  a  has  the  gamma  distribution  with  density  g(* ;(m-l ),(m-l)/aj.  The 
density  of  is  given  by 


f v  (XjT5a) 

X(m) 


nl  1  _-m(x-x)/af,  .-(x-xj/al 

(i.l[:rn-m)!ae  L1'6  J 


n-m 


for  x  >x  and  fv  (x;x,a)=0  for  x<x.  The  distribution  function  of  X,v  is 
m  X(m)  -  (m) 

given  by 


Fv  (x;x,a)  =  1-B/e  T^a ;m,n-m+l 

x(m)  ' 


8(p;m,n-m+l )  -f  r™-1  (l-r)""”  dr 

0 


where 


n 


Let  Since  and  a  are  independent,  for  -oo<z<c» 


Pt,a(ItMaiX(m)*2  a) 


*|  PT,a(x(Tn)iT+Ma-2t)fa(t;t*al  dt 


00  M 

=  J  B(e  ezt;m,n-m+1 )g(t;m-l ,m-1 )  dt 
0 


Now  =  x  +  a  log(l/q),  so  that 


oo 

pT>a(xqlx(m)+z  a)  =  J  B(qezt;m,n-m+1 )g( t ;m-l ,m-1 )  dt 


Given  0<a<l ,  define  z  =  z  (q;m,n)  by 

a  a  ^ 


®  z  t 

J  B(qe a  ;m,n-m+l )g(t;m-l ,m-l ) dt  =  a 


(2.2) 


Then 


pTjU(xqlx(m)+za  a.)  *  a  for  -=o<t<°o  and  a  >  0  .  (2.3) 


M  iq'X("l)+Z«/2  S’  VX(m)+Zl-(a/2)  a  and  I=[iq-Xq]'  ^  Eq'  (2’3) 


p  _  (x_j^ I )  =  a  for  -oo<T<«  and  a  >  0 
T  »o  q 
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In  other  words,  T  is  a  classical  100(1 -a)%  confidence  interval  for  xq  within 
the  context  of  the  two-parameter  exponential  model.  It  is  called  the  two- 
parameter  exponential  confidence  interval  procedure.  The  procedure  is  location 
and  scale  invariant.  That  is,  if  x(] ) »• • • ’X(m)  are  replaced  by  d+bX^ ^ ,. . . ,d+bX^m 
with  b>0,  then  the  left  and  right  end-points  x^  and  xq  of  I  are  replaced  by 
d+bXq  and  d+bxq,  respectively. 

Formula  (2.3)  and  hence  the  confidence  interval  procedure  I  have  a 
generalized  Bayesian  interpretation.  To  see  this,  let  ir  denote  the  generalized 
prior  density  on  (-oo,oo)x(0,«)  defined  by  TT(x,a)  =  l/a.  The  corresponding 
posterior  density  is  defined  by 


TT  (  t  » a  |  x  i ,. . . ,  x^ ) 


for  x  <  x_  and 
—  m 

7r(x,ajx1 ,. . .  ,xm)  =  0  for  x  >  xm 

(It  is  understood  here  that  x, >...>x  .) 

I  m 
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Let  ^  and  ^  denote  random  variables  having  joint  density  fffx.alx-j .  ,x  ). 
Then  the  marginal  density  Tra(a|  x] . .  ,xm)  of  ^  is  given  by 

TTa(a|  xi  ....  ,xm) 

% 

CO 

j  tt(t ,a | x-j ....  ,xm)  dx 


00  CO 


J  da  f  tt(t ,a  I  x1 . .  ,xm)  dx 

0  —CO 

=  \  9  ^|;rn-l  ,m-l^ 

where  a  =S'J5(xi-xm)/(m-l )  and  g(*;m-l,m-1)  is  a  gamma  density  as  defined 
above.  The  conditional  density  1T^|^(Tla;xi  »•  •  •  »xm)  of  ^  given  £  =  a  is 


obtained  as 


,  1r<T.a|x1,....xm) 

00 

J  if(x,a | x-j » . . .  ,xm)  dx 


n!  -1  ,-m(xm-T)/ar  -(VT)/al 
(m-1 ) !  (n-m) !  a  e  ' 


n-m 


for  x  <  x  and  tt  i  a(t  |y;x, .  ,x_)  =  0  for  t  > 
ns  t  a  i  nn 


'va, 


i 

i 

! 
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-(x  -x)/a 

Set  £  =  e  and  note  that  £  =  xm  +  ^  log  The  conditional 

density  ir  .  ,(p|  a;x-, . .  ,xm)  of  p  given  a  =  a  is  obtained  as 
p  |  a  i  in  <v 


"^p  |  a  ( p  I y  *  X1  *  *  *  *  *  xm) 


V  % 


=  f\|a(xm+a  log  Plasx! . Xh,) 


'V'Vf 


n!  r  pm_1(l-p)n"m  for  0<P<1 


(m-1 ) 1 (n-m) 1 


and  tt  ,(p|a;x-| .  ,xm)=0  otherwise.  Consequently  the  conditional  distribution 

function  of  p  given  a  =  a  is  obtained  as 
%  ^ 

F^| a^pl a;xl ’  —  *xm^  =  B(P;m>n-m+1) 


For  0<q<l  set 


Sq  =  l  +  %  log(l/q)  =  xm  +  £  log(^/q) 

Then  the  conditional  distribution  function  of  x.  given  a  =  a  is  obtained  as 

q 

Fx  la(xla;xl . xi») 

/  (x-x  )/a  \ 

=  We  iaix’ . v 

/  (x“xm)/a  \ 

=  B (q  e  ;m,n-m+l  ) 


I 
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B 


Consequently  the  distribution  function  of  is  given  by 


Therefore 


F  (x|x, - - xm) 


Q 

( 


Fx  | a( x |  a;xi , -  - .  ,xm)7T  (a  |  x1 , . . .  ,xm)  da 


^q 


/  ( X-  X  ) / 3  \  a  /a  \ 

=  J  B|q  e  ;m,n-nH-l  j  g  ^|-;m-l  ,m-lj  da 


00  /  (x-xjt/a 


n  V 


q  e  ;m,n-m+l )g(t;m-l ,m-l )  dt 


Fxn(xm+zaaixl,",,xm^  s  “ 

where  z^  is  defined  by  Eq.  (2.2).  This  yields  the  generalized  Bayesian 

interpretation  of  Eq.  (2.3). 

2.1.1  Modifications  to  Handle  Grouped  Data 

In  many  applications  X-j,...,X  are  rounded  up  or  down  or  grouped  to 

yield  a  small  to  moderate  number  of  distinct  values.  This  rounding  or 

grouping  can  have  an  adverse  effect  on  the  confidence  interval  procedure 

described  above  unless  the  procedure  is  appropriately  modified. 

Specifically  let  k>2  denote  a  positive  integer  and  let 

-<«=dg<d^<. .  .<d^  i<d^=<».  Fo:  1  <j < k  let  Nj  denote  the  number  of  sample 

values  in  the  interval  (d.  -,,d.).  Then 

J-l  J 
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P 

t,a 


(Ni 


=n 


1 


.Hk-"k) 


=  C(n1 . ,nk)  ^ [F(dj-,T,a)  -  F(dJ._1  ;x ,a ) ]  3  , 

where  F(-«>;T,a)=0,  F(°°;x,a)=l ,  n-|,...,nk  are  nonnegative  integers  adding  up 
to  n  and 


C(n-| , . . .  ,nk)  ~ 


n: 


'r 


•  n. 


If  n^>0,  then 


PT,a(Nrnl’-,,,Nk=nk) 

n,  n . 

-(d-| -x)/a1  k  r-(d..rx)/a  -(dj-T)/a“l  J 


L(n 


r  -(d,-x)/a1  1  k 

,nk}  |l-e  J  e 


n . 


=  C(n, , . . .  ,n,  )  77*  e 

j-2 


k  |  -(dj_1-d1  )/a  -(dj-d-| )/a 


3  -(n-n-j )  (d^  -x)/a 


1-e 


-(d-|  -x)/a 


-."l 


for  x  <  d-j  and  a  >  0,  while  the  indicated  probability  equals  zero  for  x  ^  d-j . 

An  approximate  100  (l-a)%  confidence  interval  for  will  now  be  obtained 
by  modifying  the  generalized  Bayes  derivation  of  T  given  above.  Let  tt  again 
denote  the  generalized  prior  density  on  (-<*>, °°)x(0,°°)  defined  by  ir(x,a)=l/a. 
Suppose  l<n^<n.  Consider  the  corresponding  posterior  density  defined  by 
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tt ( t , a  1  n -j . .  ,nk) 


'  Tr(T»a)pT.a(|Vni  *Nk=nk^ 

"ZrC^(T,a)PT'a(N-|=nj ,Nk*nk)]  dr  da 


k 

a-17r 

3  J-2 

-( d . -d, )/a  -(d.-d, )/a" 

e  J  -e  J 

ni 

J  -(n-n1 )(d1-x)/a 
e 

-(dT  -x)/a" 
1-e  ' 

nl 

ff'-'k 

"-(d.  i  -d,  )/a  -(d.-d,  )/a" 
e  J  -e  J 

J  -(n-n^ )(d^-x)/a 
e 

_  ^ 

-(d,-r)/a 1 
1-e  1 

nl 

dx  da 

for  t  <  d1  and  a  >  0,  while  ir(r,a|n^ .  .n^)  =  0  otherwise.  Observe  that 

nl 

dx 


,d1  (n-n1 )(d1-x)/a 

/  e 

-(d,-x)/a" 
1-e  1 

-Loo 

- 

r~  -(n-njx  _  n, 

=  a  |  e  (1-e  T)  dx 

0 


1  n-n,-l  n,  n,;(n-n,-l)! 

-  a  f  p  1  (1-p)  1  dp  =  a  ! 

J0 


n: 


Let  £  and  ^  denote  random  variables  having  joint  density  7r(x,a  jn-j , . . .  ,nk 
Then  the  marginal  density  tt  (a|n1 .  ,nk)  is  given  by 


a* 

lL-(dj-rV'a  -Cj-HWf1 


n. 


k  ' 
7T  e 

.  tL L 


-e 


r k  r 

/  7 T  € 

JQ  J=2  L 


k  r-(dj_1-d1)/a  -{dj-d1)/arinJ 


(2.4) 


-e 


da 


The 


conditional  density  irx|  (x|a;ni ,. . . ,nk)  of  t  given  ^  =  a  is  obtained 


as 


^|^(Tla;n1....»nk) 


r 

n,  -(n-n])(d1-T)/a  I"  -( d-j ~T)/a~] 

n-j  (n-n-j-1 ) !  e  L1  "  e  J 


for  t  <.  d-j ,  and  | ^ ( t | a ; n -j  ....  ,nk)  =  0  for  t  >  d-| . 

Set  £  =  and  note  that  £  =  d-j  +  a  log  £.  The  conditional 

density  TT^| a(p| a;n1 .... ,nk)  of  £  given  ^  =  a  is  obtained  as 

Tr^U(p|a;nl,'",nk) 

=  p^l^l  +  a  l0g  Pi  aini . %) 


n'  "-Hr1  ni 

n -j 7"( n - ri -j - 1 J !  p  (]-P)  for  0  <  p  <  1 


and  it  (p|  a;n] .  ,nk)  =  0  otherwise. 

%  I'v 


The  conditional  distribution  function  of  £  given  a  =  a  is  obtained  as 
|j(p|a>0-j,...,nk)  -  B(p;n-Oj  ,n^+l ) 


where  B(p;a,3)  is  as  defined  above. 
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For  0  <  q  <  1  set 


xq  =  l  *  £  1og  O/q)  5  dl  +  ^  lo9  (j^/Q) 

Then  the  conditional  distribution  function  of  xq  given  ^  =  a  is  obtained  as 

'WXla''"l . "k> 


/  (x-d, )/a  \ 

F*u(qe  |a;ni . "k) 

/  (x-d, )/a  \ 

=  Blq  e  ;n-n1  .n^+l J 

Consequently  the  distribution  function  of  is  given  by 

Fxq(xh . nk> 


-  f  Fx  |  (x|a;n1 . ,nk)  t  (a |n] , . . .  ,nk)  da 
jq  q 1  “  K> 


“  /  (x-d, )/a  v 

=  J  B|q  e  ;n-n1  .n^l  j  ir^(a|n1 .  ,nk)  da 


where  -r  (a  |  n^ , . . .  ,nk)  is  given  explicitly  in  Eq.  (2.4). 


Given  0  <  i  <  1,  define  x^  =  x<^(n1 , . . .  ,nk)  and  x^  =  xq(n1 , . . .  ,nk)  by 


Fx  =  2  and  Fx  (xqlnr---,nk^  1  '  2 


~3 


^3 
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Also  set 

I  =  Xq  (N-| , .  • .  ,Nj^)  J 

Then  T  should  be  a  good  substitute  for  the  two-parameter  exponential  100  (l-a)% 
confidence  interval  when  the  sample  data  are  rounded  or  grouped.  (Note  that 
additional  rounding  or  grouping  may  be  required  to  guarantee  that  0  <  N-j  <  n. 
This  provides  no  problem  in  practice,  for  if  the  requirement  cannot  be  met, 
the  original  rounded  or  grouped  observations  are  identical  and  there  are  no 
reasonable  confidence  intervals  for  x^.) 

2.2  TWO- PARAMETER  WEIBULL  PROCEDURE 

Consider  the  two-parameter  Wei  bull  model 

F(x;t,B)  =  1  -  e_t' 

=  0 

where  t  >  0  and  0  >  0.  Correspondingly 

-tx0 

S(x;t,8)  •  e  zx  x  >  0  , 

=  0  x  <  0 


x  >  0 

x  <  0 


0  <  q  5  1 
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-1  /S 

In  this  model  8  is  a  shape  parameter  and  a  =  t  is  a  scale  parameter. 
Alternatively  log  a  and  1/8  are  respectively  the  location  and  scale  parameters 
for  the  distribution  of  log  X,  where  X  has  distribution  function  F(*;t,8). 

Given  2  <_  m  <_  n,  let  X^,...,X^  denote  the  upper  m  order  statistics 
based  on  a  random  sample  of  size  n  from  F(*;t,6).  The  joint  density  of 
these  random  variables  is  given  by 


<!>...  ,xm;t,3) 


The  maximum  likelihood  estimators  of  t  and  8  are  not  easily  found. 

An  approximate  generalized  Bayes  confidence  interval  procedure  for  xq 
will  now  be  obtained.  Let  tt  denote  the  generalized  prior  density  on 
(0,°°)x(0,<»)  defined  by  ir(t,3)  =  l/t8.  The  corresponding  posterior  density 
is  defined  by 

tt(  t,8|  x-j , . . . , x^ ) 
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yB(x, 


f  dB  f  tt( t , B| x,  .  ,x  )  dt 
0  0 


A  /  JJ)  \  ^  OO  ^ 

(t8)m_1  xi  j  J  g(tx®,t/t)  dt 


°°a  ,  /  rn  \  °°a  /•. 

J  ( t8)m_  (^TTxiJ  dB  J  g(tx®,t/t)  dt 


The  conditional  density  -rrt ,  ( 1 1  g;x^ , . . .  ,xm)  of  t  given  £  =  6  i  s  obtained  as 


tt( t , S { x-j . .  ,xm) 

f  -rr(  t ,  S  |  x-j , . .  -  ,xm)  dt 
J0 

g(tx®,t/t) 

00  A  A 

J  g(tx^,t/t)  dt 
0 

The  corresponding  conditional  distribution  function  of  £  given  £  =  8  is 


obtained  as 
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F  (x|x,,...,x  ) 

— q 

oo 

=  jf  Fxql^^X^  6;X1  ’  *  "  ,Xm^  ^(2 !  x-j , - . .  ,xm)  de 

Given  0  <  a  <  1  defined  and  x^  in  terms  of  X^,...,X^  according 
to  the  formulas  Fx  (  x^l  j ,. . .  ,X^)  =  a/2  and  Fx  (  xq|X^  j .  .X^j)  =  1  -  a/2. 

Also  set  I  =  [x  ,x  ].  Then  P.  „(x  il)  =  a  for  all  t  >  0  and  8  >  0.  The  proof 
~q  q  t ,  fcj  q 

of  this  result,  which  will  not  be  given  here,  depends  on  the  observation  made 
at  the  beginning  of  this  subsection  that  the  two-parameter  Weibull  model  can 
be  viewed  as  a  location-scale  model.  The  result  shows  that  T  is  a  classical 
1 00( 1 -a)%  confidence  interval  for  xq  within  the  context  of  the  model.  It  is 
called  the  two-parameter  Weibull  confidence  interval  procedure.  The  procedure 
is  invariant  under  scale  and  power  transformations.  That  is,  if  X^,...,X^ 
are  replaced  by  d  Xb-^,...,d  Xb^  respectively  with  d  >  0  and  b  >  0,  the  left 
and  right  end-points  x  and  x  of  I  are  replaced  by  d  xD  and  d  jr  respecti vely. 

i  *  n  Q 

[To  see  this,  observe  first  that 

tt( t , 0 1 d  xb,...,d  xb)  =  dB  bTT(d8t,b6’x1 ,. . .  ,xm) 

and  then  use  this  equation  to  show  that 

Fx  (d  xb|d  xb,...,d  xm)  =  Fx  (x jx] ,. . . ,x  )  .] 

The  confidence  intervals  for  x^  can  be  transformed  as  described  in  Subsection  2.1 
to  yield  100(l-a)%  confidence  intervals  for  S(x). 
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Unfortunately  the  two-parameter  Weibull  confidence  interval  procedure 
is  not  computationally  feasible.  To  obtain  an  approximate  version  of  the 
procedure  which  is  computationally  feasible,  for  v  >  0  and  t  >  0  set 
h(v,t)  =  log  g(v,t)  =  (m-1)  log  t  +  (n-m)  log  (l-e~vt)  -  mt  and  let  h'(v,t), 
etc.,  denote  differentiation  of  h  with  respect  to  t.  Then 


h'(v.t)  .stl  +  telv. 

1  e  -1 


b  (v»t)  = 


m-1  (n-m)v2evt 
.2  .2 
C  (ejt-l ) 


tT'lv.t)  =  +  IsjflvVVM)  >  o 

1  <e'*-l) 


Thus  h(v,t)  has  a  unique  maximum  at  t^  =  t(v)  which  is  the  unique  root  of 
h'(v,t)  =  0.  This  root  can  be  found  by  applying  Newton's  method  to  the 


function  h' (v,* ) . 


Consider  the  approximation 


h(v,t)  •  h(v,t  )  +  j  h"  (v,tQ)(t-t0)2 


=  log  g(v,t  )  +  j  h  "  (v,t  )(t-t  )2 


Correspondingly 


h"  (v,t  )(t-tJ2/2 
g(v,t)  i  g(v,tQ)  e  00 


This  can  be  written  as 


g(v,t)  i  C(v)N(t;t(v) ,a  (v) ) 


where 


a(v)  =  \-l/h  "  (v,t(v) ) 


C(\j)  =  g  (v,t(v)  )a(v)  V2 7  , 

2  2 
and  N( • ;p,a  )  is  the  normal  density  with  mean  y  and  variance  a  .  The  approx¬ 
imation  for  g  in  turn  yields  the  approximations 


tm6m-1^x1j  CM3)) 

I  tm  6m_1  (r\  x^  C(C(3))d6 


irt  g(t|B;x1,...,xm)  i  N  (t;t  t(v(B)),  t2  a2(v(  3))) 


A  A  Q 

where  v(b)  =  t  x^. 
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Let  $  denote  the  standard  normal  distribution  function  and  set  Q  =  1  -  4>. 
Then 


(t|6;x-| 


’V 4  i 


t  t-t(  v(B))  \ 
,  cr(v(6))  / 


so 


Fx  |e(x|B;x,,...,x  )  5  q(^!Jo2UM±<2M 

1  ”  V  o(C(S)) 


Set 


Fx  (xlxl 


’xm} 


,Xm) 

m 


dB 


> 


where  ir0  and  F  lo  are  the  approximations  to  and  F  ,Q  just  determined. 

6  13 

Given  0  <  a  <  1  define  x„  and  x  in  terms  of  X/, \ . . ,X/  x  according 

-q  q  (I)’  (m) 

to  the  formulas 


a 


and  F 


X  ' 

JLq 


Xq'X(l)’- 


a 

2 


Then  I  =  [x^x^]  determines  the  approximate  two-parameter  Wei  bull  confidence 
interval  procedure.  It  is  also  invariant  under  scale  and  power  transformations. 
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2.3  QUADRATIC  TAIL  PROCEDURE 

Let  3  <  m  <  n  and  let  x(] ) » • • • >x(m)  denote  the  upper  m  order  statistics 
based  on  a  random  sample  of  size  n  from  F.  Set  £(x)  =  -log  S(x)  =  -log  [l-F(x)]. 
Then  S(x)  =  e‘£^,  so  that 


s(x)  -  S(.X(m))  e 


■[i(x)-A(X(m))] 


x  >  X 


(m) 


In  particular,  if  xq  >_  X^,  then 


q  -  S(X(m))  e 


'^(xq)“£(X(m) )] 


and  hence 


=  l09 


S<X0.)> 

q 


In  other  words,  xq  is  the  solution  to  the  equation 

*(x)  -  A(X^mj)  =  y  ,  (2.5) 

where  y  =  log  Cs(X^m^)/q].  The  solution  to  Eq.  (2.5)  can  be  written  as 

x  =  X(m)  +  L(y)  ,  (2.6) 


where  L  depends  on  X^  as  well  as  the  (unknown)  distribution  function  F. 
By  Eq.  (2.6) 


Xq  =  X(m) 


L  log 


S(X 


(m) 


o  <  q  £  s(x(mj) 


(2.7) 
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This  suggests  estimates  xq  of  xq  having  the  form 


('-  V)  ■ 


xq  =  X!,»)  +  L('°9  '  0<"iS(lW 

It  is  natural  to  estimate  S(X^)  by  s(x(m))=  m/n-  This  leads  to 

:(109  ft)  • 


Xq  =  X(-D  *  L' 


0  <  q  <  - 
^  —  n 


(2.8) 


Consider  for  example  a  distribution  function  F  belonging  to  the  two- 
parameter  exponential  model.  Then 


F(x)  =  1  -  e'(x-T)/a  ’  x 


where  -«  <  T  <  «  and  a  >  0.  Correspondingly, 


£(x)  =  ,  x  >_  t  , 


and 


L(y)  =  ay  ,  y  >  0 


~  ,  m-1 

a  =  S-T  ^  [X0)"X(m)] 


If  a  is  estimated  by 
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and  L(y)  =  ay  by 

L(y)  =  ay  ,  y  >  0  ,  (2.9) 

Eq.  (2.8)  reduces  to  the  exponential  tail  estimator  described  in  Subsection  2.1 . 
A  natural  extension  of  Eq.  (2.9)  is  estimators  of  L  having  the  form 


L(y)  =  ay  +  |  y2 


y  >  0 


(2.10) 


~  2 

where  the  additional  term  b  y  /2  hopefully  properly  takes  into  account  small 
to  moderate  departures  of  the  tail  of  F  from  the  two-parameter  exponential 
model.  Together,  Eqs. (2.8)  and  (2.10)  yield  the  estimator. 


x  +  X/  v 
q  (m) 


+  a 


0  < 


q  <  - 
m  _  n 


(2.11) 


It  is  desirable  that  xq  be  a  nondecreasing  function  of  1/q.  This  is 

A  A 

true  for  the  estimator  given  in  Eq.  (2.11)  if  a  and  b  are  nonnegative.  Other¬ 
wise  the  estimator  can  be  modified  in  an  obvious  way  to  make  it  nondecreasing  in 

/N  A 

1/q  [set'x  =  x  for  q<qQ,  where  qQ  is  chosen  as  large  as  possible  subject  to  the 
constraint  that  L  given  byEq.  (2.10)  is  nondecreasing  in  y  for  0  ^y  <_  log  (m/nqQ)]. 

A  /V 

In  order 'to  determine  specific  choices  of  the  quantities  a  and  b 
appearing  in  the  definition  of  L,  it  will  be  assumed  that 


L(y)  =  ay  +  |  y2  ,  y  >  0 


(2.12) 


where  a  ^  0  and  -°°  <  b  <  Of  course  this  can  be  exactly  true  only  if 
b  >_  0,  so  the  following  discussion  is  "formal"  if  b  <  0. 
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By  Eqs.  (2.7)  and  (2.12) 


S(X,  *)  . 

■  X(m)  *  a  +  I 


log 


s(xM)l2 


A 

If  m  is  reasonably  large,  then  m/n  =  S(X^)  =  s(X^mj)  and  hence 


:  X(m)  ta  (“9S?) 


Set 


L  =  log  ^and  N  =  \  log  ^ 


Also  let 


9  =  a  +  Nb 


be  considered  as  an  estimator  of 


9  =  a  +  Nb 

Then  Eqs.  (2.13)  and  (2.11)  can  be  rewritten,  respectively,  as 

xn  4  X,mv  +  L9 
q  (m) 

and 

A  A 

xn  =  X,  >  +  L9 
q  (m) 

Similarly  a  confidence  interval  J"  =  [9,9]  for  9  yields  a  confidence 

T  =  [x  ,3T]  for  x„,  where  x  =  X,,  +  L9  and  x„  =  X,m ,  +  Le\ 

u-q  q  q  -q  (m)  -  q  (m) 


(2.13) 


(2.14) 


(2.15) 

interval 
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It  is  natural  to  consider  estimators  9  of  e  which  are  linear  combina¬ 


tions  of  X 


(k) 


x(m)’  1  <  k  <  m,  with  constant  coefficients.  Such  an 


estimator  can  be  written  in  the  form 


m-1 

=  I  k  u>k[X(k)  -  X(k+1)] 


(2.16) 


It  will  be  shown  in  Appendix  I  that  the  expected  value  of  such  an  estimator 
9  is  given  by 


-  /m-1  \  /m-1  \ 

E  9  =  ^  »k)«  ♦  (I  uk  »kjb  ,  (2.17) 


where 


m-1 


wk  •  I  j 

j=k  J 


Given  an  integer  J  such  that  2  £  J  £  m,  set 


SJ  =  J-l  |  ^X(k)  '  X(J)^  =  J-l  |  k^X(k)  '  X(k+1)^ 


By  Eq.  (2.17) 


It  is  easily  verified  that 


Consequently 


J-T  ^  pk  ’  1  +  UJ 


ESj=a+(l  +  uj)b 


(2.18) 


It  is  well  known  that 


'£%  J-'»9  CM-  y 


where  y  is  Euler's  constant.  Thus  for  large  values  of  m  and  J  <_  m 


.  m-1 

yJ  »  1o9  jTf  ; 


so  Sj  is  an  approximately  unbiased  estimator  of  9  if 

1  +  log  =  N 


or,  equivalently,  if 


J  -  1  ;  (m-1 )e 


(2.19) 
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Suppose  from  now  on  that  Eq.  (2.19)  holds.  It  is  shown  in  Appendix  I  that 


Var(Sj)  =  [a  +  (1+N)b]2  +  (1+Xj)b2 


[2.2 0) 


where 


X 


J 


i  1 


m-1 


Equivalently 


Var(Sj)  =  3^-  [(0+b)2  +  (l+XjJb2] 


(2.21) 


Suppose  now  that  (Sj  -  0)/SD(5\)  has  approximately  the  standard  normal 

distribution.  Given  0  <  a  <  1  choose  z  /0  such  that 

x/  2 


1 


-x2/2 


dx  =  | 


Then 


P(-ZV2  S0(SJ>  iS0  -  "  SD(Sj))  :  1  -  .  .  (2.22) 

The  inequality  inside  the  probability  in  Eq.  (2.22)  can  be  written  as 

Var(5J> 


(2.23) 
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Set  Y  =  Z^j2/Jj^\  .  By  Eq.  (2.21),  Eq.  (2.23)  can  be  rewritten  as 

(Sj  -9)2  <  y2[(0+b)  +  (1+Xj)b2] 

This  inequality  can  be  solved  to  yield  the  interval  0“  <_  9  £  9+  where 


S  i  +  Y2  b  +Y>/(Si+b)2  +  b2(l-Y2)(l+X1) 


(2.24) 


In  order  to  get  calculable  estimates  it  is  necessary  to  replace  b  in 
Eq.  (2.240  by  an  estimator  b.  By  Eq.  (2.18) 


b  = 


(2.25) 


is  an  unbiased  estimator  of  b.  This  leads  to  the  confidence  interval  [9_, 6] 
for  9,  where  9.  and  9  are  given  by 


(a  a }  s 


sj  +  Y  b 


±y^(S  i+b)2  +  b^(l-Y^)(l+\,) 


(2.26) 


Let  I  =  ]  denote  the  corresponding  confidence  interval  for  determined 

bv  x  =  X,  \  +  L  e  and  x  =  X,  ,  +  L  9.  This  is  called  the  quadratic  tail 
-q  (m)  -  q  (m)  J - 

confidence  interval  procedure.  It  is  location-  and  scale-invariant.  That  is, 
if  X^  j ,. . .  ,X^  are  replaced  by  x  +  d  X^,...,x  +  d  X^,  where  d  >  0,  then 
x  and  x  are  replaced  by  x  +  d  x  and  t  +  d  7  .  [A  modification  to  this 

-q  q  -q  q 

Drocedure  in  which  Xj  is  replaced  by  zero  in  Eq.  (2.26)  will  be  discussed  in 
Subsection  2.4.2.] 
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2.4  MONTE  CARLO  EXPERIMENT 

A  Monte  Carlo  experiment  was  designed  and  run  to  compare  the  performance 
of  the  confidence  interval  procedures  described  in  Subsections  2.1  through  2.3. 
2.4.1  Experimental  Design 

The  experimental  design  is  a  modification  of  the  one  used  in  Breiman, 
Stone  and  Gins  [2],  Twenty  underlying  distribution  functions  were  considered. 
They  are  conveniently  defined  in  terms  of  four  groups,  each  having  five 

a.  u  x.  l. 

distribution  functions.  Let  F . ^  be  the  distribution  function  in  the  iL 

group.  They  are  determined  by  means  of  a  common  prescription.  Given  i,  a 

distribution  function  F .  of  a  positive  random  variable  X.  is  chosen  as  are 

five  positive  constants  b,  =  b..,  1  £  j  1  5.  Then  F. .  is  defined  to  be  the 

J  J  1/b, 

distribution  function  of  the  random  variable  J.  The  distribution  function 
Fi  and  five  constants  in  the  various  groups  are  determined  as  follows: 

1)  (Weibull)  F^  is  the  standard  exponential  distribution  function 

defined  by  F,(x)  =  1  -  e”x  for  x  >  0  and  F,  (x)  =  0  for  x  £  0,  while  b,  =  .5, 

b  1 

b0  =  .75,  b.,  =  1 ,  b.  =  1.5  and  bc  =  2.  [Note  that  F..(x)  =  1  -  e~X0  for 
c  o  4  o  1J 

x  >  0  and  1  £  j  £  5.] 

2)  (Mixed  Weibull)  F^  is  defined  by  F^(x)  =  [F^(x)  +  F^(x/5)]/2, 

while  b-j  =  .6,  b^  =  .84,  b^  =  1.04,  =  1.38  and  b^  =  1.65. 

3)  (Lognormal)  F^  is  defined  by  F^(x)  =  $(1og  x)  for  x  >  0  and 

F^(x)  =  0  for  x  £  0,  where  4>  is  the  standard  normal  distribution  function, 
while  b-j  =  .81,  b2  =  1.37,  b3  *  2.11,  b4  =  4.56  and  b£  =  10.81. 

4)  (Mixed  lognormal)  F^  is  defined  by  F^(x)  =  [F3(x)  +  F^ ( x/5 ) ]/2 , 

while  b-|  =  .88,  b^  =  1.41,  b^  =  2.01,  b4  =  3.52  and  b^  =  5.60. 
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To  explain  the  choice  of  the  constants  b..  . ,  let  the  tail  heaviness 
of  a  distribution  function  F  having  upper  .1-quantile  x  -j  be  defined  as 
"*"(x  j )/ll‘  (x  -j )]  ,  where  2(x)  =  -log[l-F(x)].  As  pointed  out  in  Reference  2, 
the  tail  heaviness  is  a  reasonable  measure  of  the  departure  of  the  upper  tail 
of  F  from  exponential  form.  It  equals  zero  for  distribution  functions 
belonging  to  the  two-parameter  exponential  model  and,  roughly  speaking,  is 
positive  for  distribution  functions  having  heavier  (i.e.,  more  slowly 
decreasing)  upper  tails  and  negative  for  distribution  functions  having 
lighter  tails.  The  heaviness  of  the  Weibull  distribution  functions 
Fn  » • • • »Fi 5  defining  the  first  group  are  respectively  .43,  .14,  0,  -.14  and 
-.22.  The  constants  b,  . ,  1  £  j  <_  5,  were  chosen  so  that  the  five  Weibull 

*  J 

distributions  provide  realistic  approximations  to  a  variety  of  data  that  have 

arisen  in  a  number  of  air  pollution  studies.  The  constants  b..,  2  <  i  <4 

^  J 

and  1  <  j  <5,  were  chosen  so  that  F. .  has  the  same  tail  heaviness  as 

1  J 

V 

The  sample  size  n  took  on  the  values  100,  200,  and  400.  Given  n 
and  the  underlying  distribution  function  F  =  F^ ,  purported  50%  and  90" 
confidence  interval  procedures  I  for  x^n  were  compared  with  respect  to 
PpUq  <  I)»  Pp(xq  >  I),  Pp(xq  t  I)  and  Ep(|I|/x  ).  These  quantities  were 
estimated  by  averaging  over  600  replications. 

2.4.2  Results 

Two-parameter  exponential,  two-parameter  Weibull  and  quadratic  tail 
confidence  interval  procedures  were  compared  for  various  values  of  m  selected 
more  or  less  by  trial  and  error  (only  some  of  which  will  be  presented).  Since 
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the  results  for  the  purported  50%  confidence  interval  procedures  were 
qualitatively  so  similar  to  those  for  the  corresponding  90%  confidence 
interval  procedures,  only  the  results  for  the  latter  procedures  will  be 
discussed. 

Although  the  performance  of  the  various  procedures  depends  somewhat 
on  which  of  the  four  families  the  underlying  distribution  function  F. .  belongs 
to  (indexed  by  i),  it  mainly  depends  on  the  tail  heaviness  of  the  distribution 
function,  indexed  by  j.  For  this  reason  and  for  simplicity  the  results  are 
presented  after  being  aggregated  by  averaging  over  the  four  distribution 
functions  having  each  given  tail  heaviness  (i.e.,  by  averaging  over  the  four 
values  of  i  for  each  j).  An  overall  aggregation,  obtained  by  averaging  the 
results  over  all  twenty  underlying  distribution  functions,  is  also  presented. 

See  Tables  1  through  3. 

Table  1  summarizes  the  results  for  n  =  100.  In  column  one  the  confi¬ 
dence  procedure  being  used  is  described.  In  column  two  the  tail  heaviness  is 
indicated  by  noting  the  value  of  the  shape  parameter  b^  which  yields  a 
Wei  ball  distribution  having  the  given  tail  heaviness.  Thus  the  shape  parameters 
.5,  .75,  1,  1.5  and  2,  respectively,  correspond  to  the  values  of  .43,  .14,  0, 
-.14  and  -.22  for  tail  heaviness.  The  overall  average  is  indicated  by  AVG. 
Column  three  (%  L  is  short  for  %  Left)  shows  the  indicated  average  of 

pF  (xq  <  T)  x  100% 

rounded  off  to  the  nearest  integer  (for  simplicity  and  to  provide  a  realistic 
indication  of  accuracy).  Similarly  column  four  (%  R  is  short  for  %  Right) 
shows  the  indicated  average  of 


Pp  (xq  >  I)  x  100% 
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Purported  90%  Confidence  Intervals  for  x^n 
TABLE  1.  SUMMARY  STATISTICS  FOR  n  =  100 


Procedure  Shape  %L  %R  %  Length 


W(1 5) 


Q  ( 40 ) 


E(15) 


E(10) 


.50 

6 

7 

13 

168 

.75 

6 

7 

13 

94 

1 .00 

6 

7 

13 

65 

1.50 

6 

7 

13 

40 

2.00 

6 

7 

13 

28 

AVG 

6 

7 

13 

79 

.50 

4 

9 

12 

125 

.75 

3 

7 

10 

96 

1 .00 

2 

8 

10 

75 

1.50 

3 

9 

11 

51 

2.00 

3 

9 

12 

38 

AVG 

3 

8 

11 

77 

.50 

7 

19 

27 

77 

.75 

6 

9 

14 

67 

1.00 

5 

5 

10 

57 

1.50 

5 

2 

8 

43 

2.00 

6 

2 

7 

33 

AVG 

6 

7 

13 

55 

.50 

10 

12 

22 

101 

.75 

7 

7 

14 

80 

1 .00 

5 

5 

10 

65 

1.50 

4 

3 

7 

45 

2.00 

4 

3 

7 

35 

AVG 

6 

6 

12 

66 

.50 

5 

13 

18 

95 

.75 

3 

5 

8 

82 

1 .00 

3 

2 

5 

71 

1.50 

2 

1 

4 

53 

2.00 

2 

1 

3 

41 

AVG 

3 

5 

8 

68 

E ( 1 5 ;  95%) 
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Purported  90 %  Confidence  Intervals  for  x^n 
TABLE  2.  SUMMARY  STATISTICS  FOR  n  =  200 


Procedure  Shape  %l  % R  %  Length 


* 


W(20) 


Q(  60) 


E(15) 


E(10) 


.50 

5 

7 

12 

125 

.75 

5 

7 

12 

74 

1 .00 

5 

7 

12 

52 

1 .50 

5 

7 

12 

32 

2.00 

5 

7 

12 

23 

AVG 

5 

7 

12 

61 

.50 

4 

7 

12 

104 

.75 

3 

6 

9 

80 

1 .00 

3 

7 

10 

62 

1 .50 

3 

9 

12 

42 

2.00 

3 

10 

13 

30 

AVG 

3 

8 

11 

64 

.50 

7 

16 

23 

73 

.75 

8 

8 

14 

61 

1.00 

5 

5 

11 

51 

1.50 

5 

3 

8 

37 

2.00 

6 

2 

8 

28 

AVG 

6 

7 

13 

50 

.50 

8 

11 

20 

94 

.75 

6 

6 

12 

72 

1.00 

6 

4 

10 

58 

1 .50 

5 

3 

8 

40 

2.00 

4 

3 

7 

30 

AVG 

6 

5 

11 

59 

.50 

5 

11 

15 

90 

.75 

3 

4 

8 

75 

1 .00 

3 

3 

5 

63 

1.50 

2 

1 

4 

45 

2.00 

2 

1 

3 

34 

AVG 

3 

4 

7 

61 

E(15;  95%) 
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Purported  90%  Confidence  Intervals  for  Xyn 
TABLE  3.  SUMMARY  STATISTICS  FOR  n  =  400 


Procedure  Shape  %L  %R  %  Length 


W(25) 


Q(80) 


E(15) 


E(10) 


.50 

4 

8 

12 

99 

.75 

4 

8 

12 

60 

1 .00 

4 

8 

12 

43 

1 .50 

4 

8 

12 

26 

2.00 

4 

8 

12 

19 

AVG 

4 

8 

12 

49 

.50 

3 

9 

12 

90 

.75 

3 

7 

10 

69 

1 .00 

3 

8 

11 

54 

1.50 

3 

9 

12 

36 

2.00 

2 

11 

13 

26 

AVG 

3 

9 

12 

55 

.50 

6 

16 

22 

67 

.75 

5 

8 

13 

55 

1.00 

4 

6 

10 

45 

1.50 

4 

4 

8 

32 

2.00 

4 

3 

8 

24 

AVG 

5 

7 

12 

44 

.50 

8 

11 

19 

86 

.75 

5 

7 

13 

65 

1.00 

5 

6 

10 

52 

1 .50 

4 

4 

8 

35 

2.00 

4 

3 

7 

26 

AVG 

5 

6 

11 

53 

.50 

4 

10 

14 

83 

75 

3 

5 

8 

68 

1.00 

2 

3 

6 

55 

1 .50 

2 

2 

4 

39 

2.00 

2 

1 

3 

29 

AVG 

3 

4 

7 

55 

E ( 1 5 ;  95%) 


I 
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and  column  five  shows  the  indicated  average  of 

PF  (x  l  T)  x  100% 
ij  q 

the  numbers  again  being  rounded  off  to  the  nearest  integer.  Finally,  column 
five  shows  the  indicated  average  of 

Ef  ( i I | / x  )  x  1 00% 

ij  4 

rounded  off  to  the  nearest  integer. 

Let  E(m;100(1-ct)%)  denote  the  (purported)  100(l~a)%  two-parameter 
exponential  confidence  interval  procedure  based  on  X^,...,X^  and  let 
w(m;100(l-a)%)  denote  the  analogous  Weibull  procedure.  It  was  discovered 
empirically  that  the  purported  100(l-a)%  quadratic  tail  confidence  interval 
procedure  given  by  Eq.  (2.26)  yields  coverage  percentages  typically  greater  than 
100(l-a)%  and  hence  to  unnecessarily  long  intervals.  To  correct  this  defect 
and  to  simplify  the  resulting  procedure  a  modified  quadratic  tail  procedure 
was  employed  in  which  Aj  is  replaced  by  zero  in  Eq.  (2.26).  This  procedure  is 
denoted  by  Q(m;100(l -a)%) .  Set  E(m)  =  E(m;90%),  W(m)  =  W(m;90%)  and 
Q(m)  =  Q(m;90%). 

The  results  for  W(15)  in  the  columns  of  Table  1  headed  %  L,  %  R  and 
%  are  identical  in  the  various  rows  because  of  the  power  invariance  of  the 
Weibull  procedure.  The  average  coverage  percentage  of  W(15)  is  87%.  This 
suggests  replacing  W(  1 5 )  =  W(15;90%)  by  say  W ( T 5 *, 9 2% )  in  order  to  obtain 
average  coverage  percentages  of  90%.  The  modification  would  cause  a  small 
increase  in  the  average  length  of  the  confidence  interval. 
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The  average  coverage  percentage  of  Q(40)  is  very  close  to  90%,  but 

the  average  percentage  of  time  the  true  value  lies  to  the  right  of  the 

interval  is  8%,  which  is  significantly  larger  than  the  desired  value  of  5%. 

This  suggests  that  a  better  modification  to  Eq.  (2.26)  than  replacing  X,  by 

y 

zero  might  be  to  keep  Xj  in  Eq.  (2.26)  but  adjust  Y  separately  for  8_  and  W  so 
that  the  average  percentage  of  time  the  true  value  lies  to  the  left  of  the 
interval  and  to  the  right  of  the  interval  both  equal  5%.  This  modification 
would  undoubtedly  cause  some  increase  in  the  average  length  of  the  confidence 
since  the  right  end-point  of  I  is  more  sensitive  to  changes  in  the  confidence 
level  than  the  left  end-point. 

The  suggested  modifications  to  W(15)  and  Q( 40)  would  presumably  lead 
to  procedures  having  similar  behavior.  Since  W(15)  requires  substantially 
more  computations  to  implement,  Q( 40)  appears  to  be  the  preferred  procedure. 

The  average  coverage  percentage  of  E ( 1 5 )  is  the  same  as  that  of 
W(15),  namely  87%,  so  E ( 1 5 ; 92% )  should  yield  average  coverage  probabilities 
of  very  close  to  90%.  A  more  serious  defect  is  that  for  the  heaviest  tailed 
distribution  functions,  the  true  value  lies  to  the  right  of  the  E ( 1 5 )  con¬ 
fidence  interval  19%  of  the  time.  On  the  other  hand,  the  average  length  of 
E ( 1 5 )  is  substantially  less  than  that  of  Q(40) .  This  suggests  modifying 
E ( 1 5 )  to  improve  its  average  coverage  percentage  for  the  heaviest  tailed 
distributions  at  the  expense  of  increased  average  length.  The  results  for 
two  such  modifications,  E( 1 0)  and  E(15;95%)  are  shown  in  Table  1.  Clearly 
E(15;95%)  is  the  better  of  these  two  procedures.  It  is  also  clear  that  still 
better  modifications  to  E ( 1 5 )  could  be  obtained  by  keeping  M  =  15,  keeping 
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the  left  end-point  of  the  interval  more  or  less  unchanged,  and  increasing 

A 

the  right  end-point  [by  setting  xq  =  +  z  a  for  some  z  >  z  975]- 

A  similar  analysis  can  be  made  of  Table  2  and  Table  3  for  n  =  200 
and  n  =  400,  respectively.  The  details  are  left  to  the  reader. 

2.5  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  STUDY 

The  results  of  the  Monte  Carlo  experiment  clearly  indicate  that  the 
task  of  obtaining  robust  confidence  intervals  for  the  extreme  quantile  x^n 
is  feasible  and  that. the  two-parameter  exponential  and  quadratic  tail 
procedures  are  promising  and  deserve  further  study.  But  no  definitive 
statement  can  yet  be  made  that  any  particular  procedure  is  best. 

One  attractive  procedure  that  has  not  been  tried  out  is  to  1)  use 
the  upper  (say)  [n/2]  order  statistics  adaptively  to  choose  a  positive 
number  a  such  that  the  empirical  distribution  of  -  x([n/2]) ’ *  *  * ’ 

X<([n/2]-l)  *  Xf[n/2])  is  "close"  t0  be'in9  exponential;  2)  apply  a  (possibly 
modified  version  of)  E(m)  to  the  transformed  data  ^ ,. . .  ,X^  obtaining 
an  interval  [x^.x^];  and  finally,  3)  apply  the  inverse  transformation  to 
obtain  the  confidence  interval  [)<^Xa,x^a]. 

A  similar  procedure  for  obtaining  point  estimators  of  xq  was  suggested 
in  Reference  2.  Surprisingly,  when  itwas  tried  out,  the  optimal  value  of  m  in  the 
sense  of  mean  squared  error  turned  out  to  be  m  =  n/2.  A  smaller  value  of 
m  is  probably  "best"  for  the  confidence  interval  problem.  Indeed,  it  has 
gradually  become  clear  that  the  confidence  interval  problem  differs  from  the 
point  estimation  problem  in  one  important  respect  that  is  not  readily 
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apparent— namely  that  controlling  bias  is  much  more  important  for  confidence 

interval  procedures  than  for  point  estimators  with  squared  error  loss. 

To  see  why  this  is  so  in  the  simplest  possible  setting,  suppose  that 

an  estimator  0  of  0  is  normally  distributed  with  mean  0+6  and  known  variance 
2 

a  ,  where  the  bias  6  satisfies  |0|  <_  b  for  some  known  number  b.  Then  the 

a  2  2 

maximum  possible  mean  square  error  of  0  is  a  +  b  .  Let  0  <  ct  <  1  and  let 
z^2  be  defined  so  that 


oc 

L 


-xz/2 


a/2 


V2tT 


a 

dx  =  2 


Consider  the  confidence  interval  I  =  [0  +  x-j,  0  +  x2],  where  x-j  is  chosen  as 
large  as  possible  and  x2  is  chosen  as  small  as  possible  subject  to  the  con- 

/V  A 

straints  that  P(e  <  9  +  x1 )  £  a/2  and  P(0  >  0  +  x2)  £  a/2  regardless  of 
3e[-b,b].  Then  x^  =  -z  ,2  a  -  b  and  x2  =  a  +  b  so  |T|  =  2(za/2  a  +  b). 

For  simplicity  let  a  be  chosen  so  that  z a^2  =  1.  Then  |7|  =  2(a  +  b).  For 

a  numerical  example  let  a  *  1  and  b  *  .5.  With  respect  to  the  mean  square  error 

2  2  ~ 

a  +  b  ,  9  is  exactly  as  good  as  an  unbiased  estimator  having  standard  devia¬ 
tion  /1.25  =  1.12,  but  with  respect  to  the  length  2(a  +  b)  of  the  corresponding 
confidence  interval  procedure,  the  unbiased  estimator  is  much  better. 
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3.  QUADRATIC  TAIL  APPROXIMATION 
3.1  THE  QUADRATIC  TAIL  FIT 


Given  a  distribution  function  F(x),  the  general  tail  fitting  model 
for  x  >  X(mj,  where  X(m)  is  the -m  highest  order  statistics,  starts  with 
writing  1  -  F(x)  as 


1  -  F(x)  =  [1  - 


l(x)  -  t(x(m)) 


assuming  some  parametric  form  for  l(\)  -  2,(X^),  and  then  using  this  fit 
to  estimate  extreme  values.  In  general,  a  convenient  form  for  defining  a 
tail  fit  model  is  to  write 


tw  - 1  <*(„))  ■  y  •  *  iV>  (3-°> 


then  solve  to  get 


x 


=  L(y) 


(3.1) 


The  fit  is  more  easily  defined  in  terms  of  the  L(y)  function.  For  instance, 
the  exponential  tail  approximation  is  defined  by  taking 

L(y)  *  ay 


*(x)  -  A(X(m))  =  [x  -  X(m)]/a 


» 


i 


or,  equivalently, 
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Then  assuming  m  large  enough  so  that 

1  "  F(X(m)^  =  m/n  =  P 


we  get  the  model 


1  -  F(x)  =  p  e 


-Cx 


(m) 


]/a 


This  model  has  been  used  with  success  (considering  its  simplicity)  in 
previous  work  to  get  estimates  of  percentiles  high  up  in  the  tails  of 
distributions  such  as  lognormal s,  Wei  bull's  and  mixtures  of  these.  The 
distributions  whose  tails  we  are  attempting  to  approximate  range  over  the 
exponential  class--that  is,  distributions  such  that  the  maximum  of  n  readings 

has,  asymptotically,  the  first  extreme  value  distribution.  However,  the 
asymptotics  generally  require  a  sample  size  much  larger  than  usually  available. 

For  instance,  our  interest  in  the  problem  arose  in  air  pollution  where  the 
interest  is  in  estimating  the  expected  maximum  or  second  highest  maximum  of 
365  readings. 

If  one  has  n  observations  as  data,  say  X^,...,X  ,  then  attempting  to 
estimate  the  expected  maximum  of  N  readings,  N»n,  requires  the  assumption  of 
a  parametric  model  that  is  valid  far  beyond  the  range  of  the  data.  The  range 
of  interest  to  us  is  the  range  in  which  1  -  F(x)  =  a/n,  0  <  a  <  1.  This 
range  is  of  interest  for  two  reasons.  First,  given  n  observations  X.j,...,X  , 
the  above  range  is  the  most  extreme  range  that  is,  in  some  sense,  "within 
reach"  of  the  data.  Second,  many  practical  questions  are  of  the  form.  Given 
these  n  observations,  how  can  the  expected  max  of  n  observations  be  estimated 
from  them? 
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"or  example,  tne  median  of  the  max  of  n  observations  is  penned  as 
the  solution  of 


1  -  F(x)  =  log  2/n 

Furthermore,  much  of  the  classical  extreme  value  theory  is  based  on  the 
x  value  satisfying 


1  -  F( x )  =  1/n 

Then,  even  to  be  able  to  apply  the  classical  theory,  this  point  must  be 
estimated  from  the  data. 

For  distributions  wnose  maximums  are  in  tne  domain  of  attraction  of 
the  first  extreme  value  distribution,  it  is  known  that 

lim  P(X>x+yjX>y)  -  e'x/a  ,  x  >  0 

y^co 

This  result  gives  some  theoretical  justification  for  the  exponential  tail 
approximation.  But,  as  we  mentioned  above,  a  substantial  sample  size  is 
required  for  the  tail  of  the  observed  observations  to  closely  follow  a 
conditional  exponential  distribution.  Further  discussion  of  this  issue 
will  be  given  when  the  concept  of  tail  heaviness  is  introduced  (Subsection  3.3). 

In  order  for  the  tail  exponential  to  be  a  reasonable  approximation  the 
tail  of  the  distribution  cannot  be  either  too  "heavy"  or  too  "light".  For 
example,  a  useful  intuitive  notion  of  tne  range  of  usefulness  is  that  the 
approximation  works  reasonably  for  Wei  bull  distributions 
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for  t  in  the  range  [1/2,2].  Put  another  way,  the  "curvature"  of  the  tail 
[corresponding  in  some  way  to  t'tx)]  cannot  be  too  far  from  constant  over  the 
range  for  wnich  the  tail  approximation  is  to  be  used. 

In  this  project,  a  search  has  been  made  for  models  that  will  give  a 
second  order  approximation  to  the  tail  shape,  where  the  exponential  tail  fit  is 
considered  the  first  order  approximation.  A  useful  model  is  a  quadratic  tail 
fit  defined  as  follows:  As  before,  let 

/(x)  -  z(X(m))  =  y  ,  x  >  X(m) 

x  -  x(m)  =  L(y)  . 

3ut  now  take 

L(y)  =  ay  +  |  y2  ,  y  >  0  .  (3.2) 

Our  work  with  this  model  has  shown  that  in  certain  areas  it  produces 
significant  results  in  tne  tail  estimation  problem.  The  following  subsections 
give  a  discussion  of:  1)  how  the  distribution  of  properties  of  the  model 
are  rerouted;  2!  the  concept  of  "tail  heaviness"  and  estimates  of  the  tail 
•  :  tes ;  3)  the  use  of  the  model  to  estimate  the  expected  maximum  and 

auantiles  in  the  range  specified  above;  and  A)  the  use  of  the  model  to  derive 
.  :r ‘i sence  intervals  for  parameters  such  as  the  extreme  quantiles  and  the 


maximum. 


■opppuqv 
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3.2  DISTRIBUTIONAL  PROPERTIES  OF  QUADRATIC  TAIL  ESTIMATES 
AND  AN  APPROXIMATION 

For  any  tail  fitting  method  defined  by  Eqs.  (3.0)  and  (3.1),  the 
distributional  properties  are  simplified  by  the  following  observation: 
Proposition  1.  Suppose  Eqs.  (3.0)  and  (3.1)  hold  exactly.  Let 
1  •••  -x(m)  be  the  m  highest  order  statistics.  Let 
1  •••  i)  be  the  order  statistics  from  a  sample  of 

size  m-1  from  an  exponential  distribution.  Then  the  joint 
distribution  of 

X(k)  "  X0n)  ’’  k  = 


is  the  same  as  that  of  the  variables  L(E ^ k j  ),  k  3  1 , ... .m-1. 

Proof.  The  variables  F(X^)  =  U^  have  the  distribution  of  uniform  order 
statistics.  The  model  assumptions  lead  to 


1  -  U, 


k) 


=  [1 


-L«.0t(kj  )  - 


xCW] 


or 

-  log  (1-Uk)  ♦  log  [1-U(m)]  =  t(x(k))  -  i(X(m))  . 

But  for  k  _>  m,  the  lefthand  side  of  above  has  the  joint  distribution  of 
Ep  j ,. . .  ,E^m_1  y  Hence,  writing, 

E(k)  =  ;'(X(k))  ‘  c(X(m)} 


» 
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and  using  Eq.  (3.1)  gives 

V)  '  X(m)  =  L(E(k)} 

The  derivation  of  the  distributional  properties  of  the  various 
statistics  based  on  the  quadratic  model  is  gotten  by  using  the  fact  that, 
in  consequence  of  proposition  1,  the  joint  distribution  of 

X(k)  "  X(m)  ’  k  =  1 ’ ' ‘ ’ ,m~l 


is  equal  to  that  of 

3  E(k)  +  !  CE(k)]2  ’  k  =  •  (3-3) 


Suppose,  now,  for  example,  that  we  want  to  estimate  the  quantile 
defined  by 


* 


x 


1/n 


Wri  te 


-  F(xl/n)  =  [1  -  FCX(m>1Ie 


l(xl/n)  '  *(x(m); 


so  then,  using  1  -  F(X^)  =  m/n  we  get 


log  m  =.  l(x(/n)-  »(«,.,) 


Therefore,  using  Eq.  (3.0),  Eq.  (3.1)  again  is 


2' 


xl/n  '  X(m)  =  L^log  =  log  m  +  b( 1 og  m)  / 
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resulting  in  the  estimate 

x*/n  =  *(m)  +  '°9  "  [a  +  (J2F!)  b 
Thus,  the  problem  becomes  to  find  a  good  estimate  for  the  parameter 


Similarly,  estimating  the  expected  max,  the  second  highest  maximum,  etc.,  can  al 
be  formulated  in  terms  of  finding  estimates  for  a  parameter  of  the  form 


0  =  La  +  Mb 


(3.4) 


where  L,M  are  known.  Since  is  a  linear  expression  in  a,b  then  it  is 
natural  to  look  for  estimates  of  8  that  are  linear  in  X^,  j  <_  m.  It  is 
convenient  to  write  these  estimates  as 


m-1 
-  Ik 
1 


a)k[X(k) 


-  X 


(k+1 


(3.5) 


Estimates 
of  are 
In 

derived. 


of  this  type  can  be  easily  computed  from  the  data  once  the  values 
given. 

Appendix  I,  the  mean  and  variance  of  the  estimate  of  Eq.  (3.5)  are 
The  main  results  are 


A 

E8 


(3.6) 


where 


m-1 

u.  =  I  (1/*) 

k  5.=  k 


(3.7) 
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Thus,  the  conditions  for  an  estimate  to  be  unbiased  are 


m-1  m-1 

L  =  S  ,  M  =  V  y,  cu 


1 


1 


k  k 


The  variance  of  f  is  given  as  follows:  Define 


(3.8) 


h  =  (b/a) 
1  k 


T,  -  v  l 
1=  1 


k  k  -A 


(3.9) 


.(2)  . 


m-1 


-  I  Mlc 


then 


Var(e)/ a' 


m-1 

=  I  Uk+h 


Vk+ 


hV 


+  h 


,  m-1 
2  's 


J.  u 


(2) 


(3.10) 


Note,  incidentally,  that  for 


cok  =  1  /m-1  , 

~  ,  m-1 

9  =  S=T  Y  [X(J)  ‘  X(m)] 

which  is  the  estimate  of  the  mean  in  the  exponential  tail  approximation. 

Using  Eq.  (3.10),  minimum  variance  unbiased  estimators  can  be  derived  as 
well  as  estimators  that  minimize  the  least  squares  loss  for  given  values  of 
h  =  a/b.  The  details  will  be  discussed  in  Appendix  II. 
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However,  minimizing  Eq.  (3.10)  even  with  a  small  approximation  as  done 
in  Appendix  II  [the  last  term  of  Eq.  (3.10)  is  discarded]  leads  to  calculations 

which  necessitate  a  small  computer  or  hand  programmable  calculator.  For 
this  reason  we  looked  at  an  approximation  which  is  applicable  not  only  to 
quadratic  L(y),  but  general  L(y). 

Write 

~  m-1 

e  ■  2k  U|i[L(Ek)  -  L(Ekt,)]  . 

Using  a  first  order  Taylor  expans ton 

L(E(k))  '  L(E(k+l))  =  [E(k)  ‘  E(k+1)]L  (E(k)  +  9^E(k)  "  E(k+1)D 
If  L'(y)  is  slowly  changing,  then  the  approximation 

L'(E(k)  *  0(Ek  -  W)  =  L'(E(E(k)>)  -  L'(Mk) 

seems  reasonable.  Now 

E(k)  "  E(k+1)  _  ^k^  ‘  k  “  l,...,m-l 
where  Z-j,...,Z  ^  are  independent  unit  exponentials.  Thus 

~  m-1 

9  ~  1  L'  (pk) 

/v  m-1  9  rj 

Var(e)  =  I  ^[L'(uk)]^ 


Then  the  variance  of  9  is 
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The  expectation  is 


-  m-1 

E(0)  =  1  o)k  l1  (uk) 

Suppose  that  L(y)  =  L(y,6j  depends  on  some  multidimensional  parameter  Q, 
and  suppose  that  8  is  an  estimate  of9(£).  Then  0  is  approximately  unbiased 
if 


9(e)  =  I  wk  L‘(uk,8) 

To  get  the  minimum  LSE  estimate  at  3  =  Sq,  look  at  the  square  error 

LSE  =  Var(0)  +  B2 


where 

B  ■  9(^1  -  ShL'<>vV  ’ 

To  minimize  the  LSE,  take  its  derivative  with  respect  to  getting 


-  B  L'  (yk,B0)  =  0 


so 


“k-  B/L'(vV 
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Not  unexpectedly, i f  L'(y)  is  increasing,  the  coefficients  decrease  and 
conversely.  For  L(y)  =  ay,  the  coefficients  are  constant,  giving  the 
exponential  tail  approximation.  B  can  be  evaluated  using 

B  =  0<BO)  -  S»kL’(vV 

=  0 (^q)  “  (m-1 )  B  ; 

so 

B  =  9(£Q)/m 

If  we  want  to  get  an  unbiased  minimum  variance  estimator  of  3,, 

~0 

A 

then  we  minimize  Var(e)  subject  to  0(3^)  =2a.k  L '  ( p k , SQ ) .  This  gives 
the  same  solution  as  above. 

Something  different  can  be  done  with  the  quadratic  tail  approximation 
to  produce  estimators  that  are  unbiased  over  the  entire  range  of  a,b  for 
which  the  tail  fit  is  reasonable.  The  LSE  minimum  estimator  of  0Q  =  L  aQ  +  M  bQ 
is  given  by 

.  1  \  ”  bn  L  +  Hho 

H  aotbouk  1+houk 

and  will  give  low  LSE  only  in  a  vicinity  of  h  =  hg.  The  unbiased  requirement, 
however,  can  be  put  in  terms  of  a  and  b  separately  because  of  the  linearity 
in  a,b.  That  is,  we  can  write 
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9  =  La  +  Mb 


E0  =  a^aj^+b 

and  require  that  the  ^  satisfy  both 

L  -  l  “k 


and 


M  =  l«kUk  . 

Minimizing  the  variance  under  these  two  restrictions  at  iig  gives 


wk  = 


XQ  +  X1 

(1+Vk} 


k 

2 


where  Xg  and  \-|  are  determined  by  the  constraints. 

However,  every  estimator  we  have  simulated  which  satisfies  the  two 
above  constraints,  no  matter  how  complicated  or  simple  it  is,  has  had  small 
bias  for  all  distributions  tested  except  the  heaviest  tailed  lognormal. 

We  give  some  evidence  in  Appendix  II  that  the  approximate  solution 
gives  results  very  close  to  those  of  the  exact  (well,  almost'.)  solution. 

3.3  TAIL  HEAVINESS  ESTIMATES 

In  the  previous  report  on  work  in  this  area,  a  definition  of  tail 
heaviness  was  proposed  as  meeting  certain  reasonable  requirements  (Reference  1, 
pp.  18-20).  The  characterization  did  not  consist  of  a  single  number,  but  was  a 


i 
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local  measure  of  the  curvature  of  the  distribution.  The  tail  heaviness 
H(p)  at  p,  0  <_  p  <_  1  was  defined  as  follows:  Let 

y(x)  =  a"  (x)/[Y(x)]2  ,  x  >  0 

and  set 

H{p)  =  -Y(xp) 

X.  L. 

where  xp  is  the  pcn  quantile  of  the  distribution.  For  exponential  distribu¬ 
tions,  H(p)  =0.  A  positive  value  of  H(p)  corresponds  to  a  tail  that  is 


"heavy"  relative  to  the  exponential,  and  a  negative  H(p)  to  a  tail  lighter 
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Lognormal 


.5 

.75 

1 .0 

mm 

1.5 

1.75 

2.0 

.3 

1.17 

.60 

.31 

.14 

.02 

-.06 

-.12 

.2 

1.03 

.55 

.31 

.17 

.08 

.01 

-.04 

.1 

.68 

.36 

.20 

.11 

.04 

.00 

-.03 

p  .05 

.61 

.34 

.20 

.12 

.07 

.03 

.00 

.01 

.52 

.29 

.19 

.13 

.09 

.06 

.03 

.005 

.46 

.27 

.17 

.12 

.08 

.05 

.03 

.001 

.39 

.23 

.15 

.10 

.06 

.04 

.02 

For  the  Wei  bull's  with  b  =  .75  to  1.5,  and  for  the  lognormal s  with 
b  21  1-25  there  is  a  minimal  curvature  problem  for  p  <_  .1,  The  curvature 
problem  is  more  severe  for  the  extremely  heavy  and  light  tailed  Wei  bull ' s 
and  the  heavy  lognormal s.  Since  the  asymptotic  extreme  value  theory  depends 
on  the  condition  H(p)  -*■  0,  the  size  of  H(p)  for  these  latter  distributions 
for  p  as  small  as  .001  indicates  the  inapplicability  of  the  theory  for  fairly 
large  sample  sizes. 

Using  the  quadratic  model  we  have: 

Proposition  1.  The  tail  heaviness  at  p  =  1  -  F(X^mj  )  is  equal  to 


h  =  b/a 
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Proof:  Since  l‘  (x)  -  2,(X^)  =  y,  then 


*'(x) 


Hx) 


A  ($L\ 

dy  l  dx  ) 


dy 

37 


Us  i  ng 


x 


2 


gives 


dx 

dy 


^  =  a  +  by 


Therefore,  at  x  =  X,^  or  y  =  0, 

-v(x(m))  -  b/a  =  h 


Therefore,  the  quadratic  model  can  be  used  to  get  estimates 
tail  heaviness.  For  a  simulation  we  constructed  four  estimators, 
three  were  of  the  following  type:  Denote 


of  the 
The  first 
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Then  we  used  a  split  half  approach  and  look  for  a  linear  combination 

C  S[m/2]  +  D 

which  will  give  an  unbiased  estimate  of  b.  Using  conditions  in  Eq.  (3/8)  gave 
the  result 


c  =  -°  =  Vu[m/2] 

The  selection  of  the  denominator  posed  a  problem.  Using  an  unbiased 
estimator  for  a  gave  a  noisy  estimate  for  h.  Finally,  we  decided  to  use 
Sm  as  the  denominator.  For  the  last  estimator  we  used  an  unbiased  estimator 
at  b  whose  variance  was  a  minimum  at  h  =  .25. 

For  the  simulation,  we  used  a  sample  size  of  200  with  the  Weibulls 
and  lognormals  as  above.  Three  different  values  of  m  were  selected, 
m  =  20,  30,  60.  The  first  three  estimators,  then,  had  the  form 


_ ] _ filmzil 

“W2]  \  SB 


m  =  20,  30,  60 


For  the  last  estimator  m  =  60  was  used,  and  it  is  denoted  by  h,,,. 

oU 

Altogether,  1000  runs  were  made,  each  run  generating  the  top  60 
order  statistics  of  the  14  distributions  using  Marsaglia's  "Super-Duper"  uniform 
random  number  generator,  and  using  the  inverse  transformation  to  get  the 
order  statistics  of  the  distributions  desired.  The  results  are  tabled  below, 
giving  the  average  and  the  standard  deviations  of  the  estimates. 


=r» 
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Wei  bull 


b 

.5 

.75 

1 .0 

1  .25 

1.5 

1.75 

2.0 

h  1 

AV 

.258 

.083 

-.003 

-.053 

-.086 

-.109 

-.126 

h20 

SD 

.329 

.315 

.307 

.303 

.300 

.298 

.296 

h  1 

AV 

.314 

.105 

.003 

-.057 

-.097 

-.124 

-.145 

h30  J 

SD 

.275 

.262 

.255 

.250 

.248 

.246 

.244 

h  1 

AV 

.422 

.134 

-.008 

-.092 

-.146 

-.184 

-.212 

oO  j 

SD 

.199 

.188 

.181 

.177 

.174 

.171 

.170 

*60  j 

.427 

.134 

-.006 

-.087 

-.140 

-.177 

-.204 

|  SD 

.187 

.159 

.148 

.142 

.139 

.137 

.136 

Lognormal 

.5 

.75 

1.0 

b 

1  .25 

1.5 

1.75 

2.0 

A 

h 

[  AV 

.525 

.296 

.180 

.111 

.066 

.033 

.010 

h20  J 

|  SD 

.377 

.357 

.344 

.335 

.330 

.327 

.324 

/s 

h 

|  AV 

.592 

.331 

.198 

.119 

.068 

.031 

.004 

h30 

l  SD 

.309 

.290 

.277 

.269 

.264 

.260 

.258 

h 

|  AV 

.693 

.371 

.205 

.107 

.043 

-.002 

-.035 

h60  ' 

\  SD 

.230 

.217 

.207 

.201 

.197 

.194 

.192 

*60 

|  AV 

.745 

.388 

.214 

.113 

.049 

.004 

-.030 

.273 

.218 

.190 

.177 

.169 

.164 

.161 

Since  the  denominator  is  an  unbiased  estimate  of  a  +  b,  instead  of  a, 
the  estimates  tend  to  reflect  the  values  of  H(p)  at  values  of  p  larger  than 

/v  /\ 

m/n.  For  all  14  distributions,  the  estimates  h^g  and  hgg  have  average  values 
close  to  H(.l).  The  averages  of  hgg  are  close  to  H(.05),  and  the  averages  of 

A 

hgg  are  slightly  above  the  H(.Ol)  values. 
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Of  course,  the  smaller  m  is,  the  higher  the  SD's  of  the  estimates. 

/V  -A 

In  use,  either  n^g  or  fi^g  is  preferable  in  terms  of  variance,  if  p  =  .1 

A, 

is  in  the  tail  range  of  interest  in  the  problem.  Note  that  n^g  almost 
always  has  a  significantly  lower  variance  than  hgg,  but  this  has  to  be 
balanced  against  the  difficulty  of  computing  the  coefficients. 

If  these  estimates  of  tail  heaviness  are  used  to  detect  departure 
from  exponential ity,  then,  using  a  2SD  rule  of  thumb,  the  true  value  of 
H(.l)  would  have  to  satisfy  |H(.1)|  >_  .3  before  the  departure  could  be 
reliably  detected.  Thus,  only  the  curvature  of  the  heaviest  tailed  Weibull 
and  the  two  heaviest-tailed  lognormal s  can  usually  be  detected.  However, 
these  are  the  distributions  that  cause  the  largest  absolute  errors  in  the 
exponential  tail  estimates. 

3.4  ESTIMATES  OF  THE  EXPECTED  MAXIMUM 

As  a  test  of  the  quadratic  model  and  exponential  tail  estimators,  a 
simulation  study  was  designed  to  estimate  the  bias  and  variance  of  a  variety 
of  tail  estimators  of  the  expected  maximum  of  the  14  Weibull  and  lognormal 
distributions. 

Twelve  estimators  were  computed  and  compared,  again  using  a  sample 
size  of  200,  1000  runs,  the  Marsaglia  random  number  generator  and  inverse 
functions. 

The  twelve  estimators  were  in  4  groups:  1)  the  maximum  was 
used  as  an  estimate  for  E  X^;  2)  three  exponential  tail  estimates 
corresponding  to  m  =  20,  30,  60;  3)  four  quadratic  minimum  variance 
unbiased  tail  estimates  with  the  variance  minimized  respectively  at 


65 


h  =  .5,  .25,  0.0,  -.15  using  m  =  60;  4)  four  quadratic  minimum  squared  error 
estimates  with  the  minimization  carried  out  at  h  =  .5,  .25,  0.0,  -.15  and 
using  m  =  60. 

The  quadratic  model  was  used  to  construct  the  third  and  fourth 
groups  of  estimates.  The  derivation  is  based  on  the  following;  Write 

X/. ,  =  X/ . ,  -  X/  ,  +  X/  >  . 

(1)  (1)  (m)  (m) 

Mow  by  the  model 

X(D  "  X(m)  =  3  E(l)  +  I  E(1) 


Use  the  fact  that 


m-1 

E ( E j )  =  y  ^  =  I  1/u 


and 


E(E?n)  =  ur 


'( 1 ) ;  '  ^1'  T  yl 


(2) 


For  m  =  60,  u]  =  4.655,  qj2)  *  1.645,  so, 


EU(i)  -  X(m^)  =  4.655a  +  1  1.657b 


Therefore,  estimating  the  parameter 

:  =  4.655a  +  11.657b 

by  -  gives  tne  estimate  X^  +  -  for  E  X^. 


Denote  by  RMSE  the  root  mean  squared  error.  The  tables  below  give 


the  value  of 

Exm 

and  the 

RMSE  of 

X^j  j  as  an 

Wei  bul  1 

estimate 

of  E  X(1). 

j 

b 

.5 

.75 

1  .0 

1  .25 

1.5 

1.75 

2.0  i 

EXo) 

36.2 

10.7 

5.9 

4.1 

3.2 

2.7 

2.4  ! 

RMSE[X(1)] 

17.2 

3.2 

1.3 

.70 

.46 

.33 

.25  i 

Lognormal 


EXd) 

63.2 

14.8 

7.4 

4.9 

3.7 

3.1 

2.7 

RMSE[X(1)] 

66.9 

8.9 

3.1 

1  .6 

.98 

.68 

.  51 

The  RMSE[X^j]  was  used  as  a  benchmark,  and  the  ratio  of  the  RMSE's 
of  all  other  estimates  to  the  RMSE [X ( -j ) 3  was  computed.  The  results  are 
tabled  below. 
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RATIOS  OF  RMSE's 


Wei bul i 


Estimates 

.5 

.75 

b 

1.0 

1.25 

1.5 

1.75 

2.0 

exp  (20) 

.61 

.61 

.66 

.73 

.80* 

.86* 

.92 

exp  (30) 

.66 

.57 

.59 

.69 

.80* 

.91 

1.02 

exp  (60) 

.95 

.59 

.49 

.71 

1  .02 

1  .30 

1  .57 

UNBMIN  (h=. 5) 

.63 

.89 

1 .11 

1.33 

1  .55 

1.76 

1  .95 

UNBMIN  (h=.25) 

.65 

.83 

.98 

1.12 

1.26 

1  .40 

1.53 

UNBMIN  (h=0.0) 

.75 

.83 

.88 

.92 

.97 

1  .01 

1.05 

UNBMIN  ( h=- .15) 

1.13 

1  .11 

1.10 

1  .09 

1.09 

1  .08 

1.08 

MINE2  (h=.5) 

.46* 

1.12 

2.47 

3.69 

4.73 

5.60 

6.35 

MINE2  ( h= . 25 ) 

.64 

.56* 

1.39 

2.23 

2.96 

3.56 

4.10 

MINE2  (h=0. 0) 

.87 

.62 

.48* 

.66* 

.95 

1.2', 

1  .47 

MINE2  (h=_J5) 

.43 

.93 

.92 

.91 

.90 

.90 

.90* 

/\a  ao  ca  on  oc  an 

Minimum  .46  .56  .48  .66  .80  .86  .90 


★ 

Minimum  value. 

Lognormal 


Estimates 

.5 

.75 

1.0 

1.25 

1.5 

1.75 

2.0 

e/o  (20) 

.48 

.53 

.54 

.54 

.55 

.57 

.59 

exp  (30) 

.52 

.56 

.54 

.51 

.50 

.51 

.52 

exp  (60) 

.61 

.68 

.61 

.51 

.44* 

.41* 

.43 

UNBMIN  (h=.5) 

.42* 

.51 

.64 

.75 

.36 

.95 

1  .04 

UNBMIN  (h= .25) 

.45 

.55 

.64 

.72 

.79 

.86 

.91 

UNBMIN  ( h=0. 0) 

.63 

.69 

.74 

.77 

.80 

.82 

.84 

UNBMIN  ('"=-.  15) 

1.17 

1.14 

1.12 

1.11 

1.10 

1.10 

1.10 

MINE2  (h=.5) 

.54 

.37* 

.42 

88 

1.37 

1  .81 

2.21 

MINE2  ( h= . 25 ) 

.48 

.52 

.35* 

.4  3* 

.67 

.98 

1  .24 

‘•'I'lE'*  (6=0.0) 

.62 

.64 

.54 

.46 

.42 

.42* 

MINE2  ( n=- .15) 

.92 

.91 

.91 

.90 

.90 

.90 

.90 

'’ini  mum 

.42 

.37 

.35 

.43 

.44 

.41 

.42 

■k 

Minimum  value. 
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The  fundamental  issue  in  the  behavior  of  these  estimates  is  the 
trade-off  between  bias  and  variance.  The  two  estimators  exp  (30)  and 
exp  (20)  have  the  best  overall  performance  with  average  ratios  of: 

Weibull  Lognormal 

exp  (20)  !  .74  .54 

exp  (30)  !  .75  .52 


The  estimate  exp  (20)  is  better  at  the  extreme  values  of  heaviness,  where 
bias  is  more  of  a  factor  while  exp  (30)  performs  better  for  moderate  to 
small  values  of  H  where  its  variance  is  smaller. 

The  other  estimators  performed  as  expected  with  one  exception,  dis¬ 
cussed  below.  The  unbiased  estimators,  as  will  be  shown  in  the  percent 
bias  tables  below,  lived  up  to  their  billing  and  had  very  little  bias  for 
any  of  the  distributions  except  the  heaviest  tailed  lognormal.  However,  the 
payment  was  in  terms  of  variance.  Except  near  the  values  of  h  at  which 
their  variance  was  minimized,  their  variances  were  large  and  produced 
inflated  RMSE's.  The  minimum  squared  error  estimates  also  lived  up  to  their 
billing  by  producing  either  the  minimum  or  near  minimum  ratio  near  the  values 
of  h  at  which  they  were  optimized.  Their  difficulty  was  that  a  bias  which 
is  small  at  one  value  of  h  can  be  large  at  other  values.  The  large  bias  led 
to  large  ratios  away  from  the  value  of  h  at  which  they  were  optimized. 

A  look  at  the  percent  of  bias,  that  is. 


100 


E  X ( 1 )  -  est 
E  X(l) 


is  revealing. 
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Notice  that  the  estimators  certified  as  unbiased  by  the  quadratic 
model  do,  indeed,  have  very  small  biases.  The  best  is  the  unbiased  estimate 
whose  variance  is  minimized  at  h  =  0.  Its  average  percent  bias  (absolute 
values)  is  2.0%.  The  only  problem,  again,  is  with  the  heaviest  tailed  log¬ 
normal  where  the  bias  rises  to  10%. 

To  show  that  even  simple  estimates  satisfying 


M  =  I  cuk  uk 

will  be  relatively  unbiased,  we  selected  an  estimate  of  the  form 

X0  S( 30)  +  X1  S(15) 

The  values  of  Xg,  X-|  satisfying  the  two  unbiasedness  conditions  are 

^  =  (m-l)/u15 

xo  =  L  -  X1  • 

24 

where  y,,  =  2  1/k.  Using  this  with  values  of  L,M  adjusted  to  estimate 
15  15 

E  X||j  with  the  top  30  order  statistics,  the  estimator  gave  the  following 
bias : 


Wei  bull  .3 


Because  of  the  fact  that  the  estimators  that  minimize  the  squared  error 

are  sensitive  to  the  choice  of  h  used,  another  experiment  was  carried  out 

for  the  distributions  such  that  H(.l)  >  0.  In  this  simulation,  the  unbiased 
/\ 

estimator  h  of  h  =  a/b  was  first  computed  using  the  upper  60  order  statistics 
then,  the  minimum  squared  error  estimate  corresponding  to  the  value  max(h,0) 
was  calculated.  The  results  are  promising,  as  given  in  the  table  below  of 
the  average  ratios  of  the  RMSE  to  that  of  X^. 

AVERAGE  RATIOS  OF  RMSE/RMSE[Xn , ]  FOR  H(.l)  >0 


Wei  bull  Lognormal 


exp  (20) 
exp  (30) 

2  step  est 


Although  the  improvement  is  not  very  large,  it  is  suggestive  for  future  work. 

3.5  CURVATURE  AND  TRANSFORMED  TAILS 

One  somewhat  strange  bit  of  behavior  is  given  by  the  estimates  optimized 
at  h  -  -.15.  On  the  one  hand,  their  biases  are  consistently  small,  doing 
better  for  large  positive  H  values  than  the  unbiased  estimates  optimized  at 
positive  values  of  h.  We  conjecture  that  this  is  due  to  the  fact  that 
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minimizing  the  variance  at  a  given  value  of  h  does  not  enhance  the  general 
unbiasedness  properties.  However,  a  harder  to  explain  phenomenon  is  that 
the  minimum  squared  error  estimate  at  h  =  -.15  does  not  give  a  significant 
decrease  in  the  squared  error  for  the  Wei  bull  distributions  with  negative 

H.  In  fact,  for  the  four  Weibulls  with  negative  H,  exp  (20)  is  a  consistent 

2 

improvement  on  MINE  (-.15)  except  for  a  slight  difference  at  b  =  2.0. 

At  first,  we  thought  that  this  indicated  some  deficiency  of  the 
quadratic  model  for  light  tailed  distributions  and  that  a  different  method 
of  selecting  the  coefficients  ^  in  the  estimate 

~  m-1 

9  =  } k  [X^>  ’  V+i)] 

would  bring  the  ratio  of  MSE's  down.  We  concluded  that  the  difficulty  was 
more  fundamental.  First  we  noticed  that  even  with  the  estimators  exp  (20), 
exp  (30),  the  ratios  go  up  rapidly  as  the  Weibulls  become  light  tailed, 
i . e . ,  b  >  1 .  On  investigating  the  cause,  it  turned  out  that  the  difficulty 
was  not  really  with  the  bias  (although  that  contributed)  but  instead  with 
the  fact  that  the  standard  deviations  of  exp  (20),  exp  (30),  were,  surprisingly 
enough,  not  that  much  less  than  the  standard  deviation  of  X^.  The  results 
are  given  below. 


SD  AND  BIAS  OF  ESTIMATORS 
WEIBULL  DISTRIBUTION 


1.0 

1,25 

1.50 

1.75 

2.0 

SD 

1.38 

.70 

.46 

.33 

.25 

bias 

0.00 

0.00 

0.00 

0.00 

0.00 

SD 

.84 

.50 

.35 

.26 

.21 

bias 

0.00 

.10 

-.11 

-.11 

-.11 

SD 

.75 

.46 

.32 

.24 

.20 

bias 

0.02 

-.15 

-.18 

-.17 

-.17 

exp  (30) 


Thus,  for  instance,  even  if  exp  (30)  were  unbiased  at  b  =  2.0,  the  ratios  of 
RMSE's  would  be  .80.  This  led  to  the  conjecture  that  the  essential  difficulty 
was  with  the  linear  form  of  the  estimate;  that  at  least  with  light  tailed 
distributions,  one  could  not  construct  a  linear  combination  of  -  ^k+-|j 
that  was  not  too  biased  without  a  resulting  SD  close  to  that  of  X^. 

This  led  to  the  idea  of  seeing  how  much  an  optimum  or  nearly  optimum 
transformation  could  improve  the  RMSE.  If  X  has  a  Wei  bull  distribution  with 
parameter  a  then  Xa  has  an  exponential  distribution.  Therefore,  given  X^, 
the  differences 


,a  v  a  ,  „  „ 

£(k)  "  X(m)  ’  k  >  m 


have  exactly  the  distribution  of  exponential  order  statistics.  That  is, 
the  transformed  variables  X^  are  exactly  fit  by  an  exponential  tail.  Thus, 
the  exponential  tail  estimate 
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e  *  X 


(m) 


I  ,  m-1 

)  1  V  TV2  yCt  -I 

U1  |  m-1  -  [X(k)  '  X(m) ■* 


should  give  a  good  estimate  of  E  y  Then  take  (e)^Xa  to  be  the  estimator 
of  E  X^j.  We  assumed  that  a  was  known  for  the  Weibulls.  That  is,  for  the 
Wei  bull  raised  to  the  power  b,  we  took  a  =  b. 

For  the  lognormal s,  it  is  not  clear  what  the  optimum  transformation  is. 
We  selected  the  powers  a  based  on  a  comparison  of  the  lognormal  heaviness  to 
the  corresponding  Weibull  heaviness.  So,  as  a  guess,  we  took  the  values  of 
a  corresponding  to  b  to  be  given  by 


b  =  .5  .75  1.0 

1.25 

1.50 

1.75 

2.0 

a  =  .4  .5  .6 

.7 

.8 

.9 

1  .0 

The  results,  expressed  as  ratios  to  RMSE[X^]  are 


Ratios 


b 

.5 

.75 

1  .0 

1 .25 

1.5 

1.75 

2.0 

Weibull 

.44 

.47 

.49 

.49 

.50 

.50 

.51 

Lognormal 

.41 

.40 

.39 

.39 

.40 

.41 

.43 

The  Weibull  results  are  a  lower  bound  under  the  assumptions  that  only 
the  top  60  order  statistics  are  used  and  that  the  scale  is  unknown.  Of 
course,  there  is  the  question  of  how  much  bias  has  been  introduced  by  the 
approximation  E  X^j  =  [E  X^]  For  The  Weibulls,  the  bias  of  the 

estimate  is  3%  at  b  =  .5  and  less  than  .6%  for  the  others.  The  lognormal 

i- 

II 
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estimates  have  considerably  larger  bias  at  the  heavier  end,  rising  to  8% 
at  b  =  1.0,  17*5  at  b  =  .75  and  39%  at  b  =  .5.  The  indications  are  that 
having  smaller  values  of  a  at  the  heavier  tailed  lognormals  would  have 
given  further  reductions  in  the  ratios. 

To  get  a  comparison  of  the  best  linear  fit  based  on  the  quadratic 
model,  we  ran  the  approximate  solution  to  the  minimum  RMSE  optimized  at 
h  =  .6,  .4,  .2,  0,  -.2.  The  table  below  is  the  minimum  ratio  over  all  of 
the  estimates. 


RATIOS  FOR  "BEST"  LINEAR  FIT 
("Best"  =  best  quadratic  model  linear  fit) 


Weibull 

Lognormal 


For  all  of  the  lognormal  distributions  the  best  linear  and  the  trans¬ 
formed  estimates  are  comparable  in  terms  of  RMSE.  They  are  comparable  for 
the  Wei bul 1  for  b  <  1.0. 

This  gives  some  evidence  that  long  and  short  tailed  distributions 
necessitate  different  procedures  to  give  good  estimates  of  E  X^.  For  the 
long  tailed  distributions  a  transformation  is  not  needed  to  "uncurl"  the  tail; 
the  appropriate  linear  combination  of  order  statistics  will  do  almost  as  well 
as  the  best  transformed  estimate.  In  the  short  tailed  cases  (Weibull  with 
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b  >  1.0)  a  power  transformation  is  essential  to  significantly  reduce  the 
SO  of  the  estimates. 

3.6  CONFIDENCE  INTERVALS 

Suppose  that  the  parameter  of  concern  is  of  the  form 

X(m)  +  La  +  Mb 

If  100%  confidence  intervals  are  found  for 

9  *  a  +  (M/L)b  =  a  +  Nb 


say  U,U,  that  is 

P0(U  <  9  <  U)  =  1  -  Q  , 

then  100Q%  confidence  intervals  for  X(m)  +  La  +  Mb  are  (approximately) 
given  by 

[X(m)  +•  LU,  X(m)  +  LU] 

Start  by  finding  a  simple  unbiased  estimator  for  9.  Take  the  estimator 
to  be  of  the  form 


S 


J 


1 

J-1 


1 


[X(k)  ‘  X(J)] 


where 
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where  J  is  the  value  to  be  determined.  Since  co^  =  1/J-l,  k  £  J-l ,  and  zero 
for  k  £  J,  the  expectation  of  Sj  is,  by  Eq.  (3.8),  equal  to 


1 


J-l 


a  +  y  uVb 


An  easy  calculation  gives 


1  J_1 

CT  y  “k  ’  1  +  “J 


Hence  J  is  determined  through  the  equation 


U,  =  N  -  1 


The  well-known  approximation 


V  l/i'vlogk  +  y 
i=l 


leads  to  Uj-log  [ (m-1 )/ (J-l )]  so  that  we  get  the  equation 


J  -  1  -  ( m- 1 )  e 


-(M-1 ) 


(3.11) 


To  compute  the  variance  of  S,,  use  Eq.  (3.10)  and  the  fact  that 


1/J-l  ,  k  <  J-l 
!/'<  ,  k  >  J 
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to  get 

Var(Sj)  =  C]  a2  +  C2  ab  +  C3  b2 

where 


m-1  ? 

^  =  I  <*>£  =  i/j-i 


a  (7^7[2(j'”  *  (J',)  uj] 

=  TPTJ  (1  +  N> 


C 


3 


( J-l  )2 


J-l  ?  m-1  _ 

I  (u  +DZ  +  1  (1/k)2 

1  *  J 


(J-D2 


m-1 

Y 

T 


.  (2) 

“k 


Now 


m-1 


I  (1/k) 

J 


2 


m-1 

V 


■J2> 


m-1  m-1 

=  I  I  1/  -  2 

1  j  =  k  J 


J-l 

Y 

T 


(pk+1)2 


( j-i ) c(2+’jj)2  +  n  -  u' 


where 


=  x  1/j 


(This  last  calculation  is  carried  out  by  writing  uk  =  Uj  +  u'k  where 
uk  =  I  1/j  • )  Hence 

j  =  k 


'3  J- 


or,  assuming  J,  say,  >_  10,  then 


tV  [(2  +  n,)2  *  1]  +  -di_2  +  u(2) 


(J-l) 


C3  =  CT  C(2  -  Uj)2  -  {0-D  u<2)] 


And  since  2  +  Uj  =  N  +  1  , 


V  ar(S 


=  JJ-T) 


pyy  ja2  +  2(N+1 )ab+  (N+l)2  b2  +  b2[l  +  (J-l )  uj2)] 


Note  that  Aj  =  (J-l)  Uj2^  -  1  -  (J-l )/ (m-1 ) . 
We  can  write  the  above  as 


Var(Sj)  =  |[a  +  0+N)b]2  +  (1  +  *j)  b2  J 


(3.12) 


The  righthand  side  of  Eq.  (3.12)  involves  the  unknown  parameters  a  and  b. 


However,  we  can  write 


a  +  (1+N)b  =a+Nb+b=e+b 


i 
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We  can  get  an  estimate  of  b  by  using  an  expression  of  the  form 

b  =  C(S,  -  S  ) 
v  J  m' 

Looking  at  Eq.  (3.8),  b  will  be  an  unbiased  estimate  of  b  if 

C  =  1/uj  =  1/N-l 


The  key  to  the  remaining  part  of  the  computation  is  the  assumption 
that  Sj  has  an  approximately  normal  distribution.  Linder  this  assumption 

P0(-z  •  SD(Sj)  <  Sj  -  9  <  z  SD(Sj))  =  1  -  Q  (3 

where  SD(Sj)  is  the  standard  deviation  of  Sj  and  z  is  computed  from  the 
normal  tables,  i.e.,  for  a  unit  normal  Z, 

P(Z  >  z)  =  Q/2 

Write  the  inequality  inside  the  probability  in  Eq.  (3.13)  as 

(Sj  -e)2  <  z2  Var(Sj) 

or  putting  y  =  z/(l-J),  and  using  Eq.  (3.12) 

(Sj  -e)2  <  Y2  [(0  +  b)2  +  (1  +  Aj)  b2]  ■ 

Simplifying  leads  to 

02(1  -  y2)  -  20(Sj  -  y2b)  -  y2  b2 ( 2  +  Xj)  +  S2  <  0 
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Solving  gives  the  following  expression  for  the  roots; 

(Sj^,2b)  *  x/(SJtr2b)2  -  Sj(l-Y2)  +  (l-Y2)Y2b2(2+Aj) 

TT7  ; 

and  simplifying  the  square  root  gives 

(Sj  +  y2b)  iY>/(2j+b)2+  b2(l-Y2)(l  +  -j) 

T77 

To  get  computable  confidence  limits,  b  is  replaced  by  the  estimate  b. 
However,  this  adds  extra  variability  to  the  limits  and  tends  to  make  them 
too  large.  To  adjust  and  simplify,  we  replace  the  factor  1  +  by  1, 
arriving  at  the  final  form 

S,  +  y2  b  +  yj{ S  T+b)2  +  b2(l -y2) 

U,U  =  — - - — - 

1  -  Y 
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4.  EXPONENTIAL  TAIL  ESTIMATES  APPLIED  TO  STATIONARY  SEQUENCES 

A  simulation  experiment  was  carried  out  to  see  how  sensitive  exponential 

tail  fitting  was  to  the  presence  of  correlation.  Recall  that  if 

X(i)  >_  are  the  order  statistics,  then  the  exponential  tail 

fit  for  x  >  X-  \  has  the  form 
-  (m) 

- [x-X/  n ]/a 

1  '  F(x)  =  1  '  F(X(m))  e 
The  parameter  a  is  estimated  by 


a 


1 

m-1 


m-1 


X(m)] 


(4.G) 


If  the  exponential  tail  fit  holds  exactly,  then 

X(k)  '  X(m)  =  a  E(k)  ’  k  = 


where  >_  . . .  >.  E^m_^  have  the  distribution  of  order  statistics  in  a 
sample  of  size  m-1  from  an  exponential  distribution  with  unit  mean.  Therefore, 


E(X 


(1) 


'  X(m)} 


E  (E 


(1 


-  a  u. 


where 


m-1 

u  =  ^  1 /k  =  log  (m-1 )  +  y 

1  1 


J 
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where  y  is  Euler's  constant.  Thus 

E  X(l)  =  E  X(m)  +  a  U1 

Estimating  E  by  and  a  by  the  estimation  a  of  Eq.  (4.0)  gives  the 

exponential  tail  estimate 

X(m)  +  U1  a 

for  the  expected  maximum  of  n  observations. 

The  stationary  time  series  were  generated  as  follows:  Let 

. ™  (4-!> 


Y 


1 


=  e 


1 


where  en  are  independent  N(0,1)  variables.  Thus,  the  Y-p.^.Y^gg  form  a 
Gaussian  Markov  chain  with  auto  correlation 


e  y  y, 

n+k  k 


oU! 


They  have  mean  zero  and  variance  one.  The  actual  sequence  used  consisted 
of  the  lognormal  variables 


Xn  = 


p  .+o  •  Y 
3  3  n 


e 


,  n  =  1  , . . . ,200 
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where,  as  in  our  previous  work, 

Uj  =  -Oog  2)/2b.  3.  =VTogT/bj 

and  b,  took  the  values  .5,  .75,  1.0,  1.25,  1.5.  This  range  included  a 

J 

very  heavy  tailed  distribution  (b  =  .5)  and  ranged  up  to  the  light  tailed 
distributions  at  b  =  1.25,  1.5. 

Each  run  was  repeated  1000  times  using  the  Massaglia  "Super-Duper*1 
uniform  random  number  generator  and  the  Box-Mullen  transformation.  Since 
the  exact  value  of  E  X^j  was  difficult  to  compute  (except  for  o  =  0)  the 
average  value  of  over  the  1000  runs  was  taken  as  the  target  figure. 
Examination  of  the  SD's  of  X^  showed  that  the  error  in  using  X^  as  an 
estimate  of  EX^j  would  not  appreciably  affect  the  results. 

As  in  our  other  work,  the  SD[X^]  was  taken  as  the  benchmark.  The 
RMS  error  at  the  exponential  tail  estimate  using  m  =  30  was  computed  using 
the  set  of  c  values 

o  =  0 ,  +  .  2 ,  +  .  4 ,  +.  8  , 

and  divided  by  SD[X^j]  to  give  a  measure  of  the  improvement  in  using  the 

exponential  tail  estimate  instead  of  X,,,  as  an  estimate  of  EX,,,. 

U)n  (1) 

First  of  all,  it  is  interesting  to  note  how  EX^  [or  rather  is 

affected  by  the  correlation.  This  is  tabled  below: 


1 

1 

b 

.5 

.75 

1.0 

.8 

48.6  1 

12.1 

6.3 

.4 

60.9 

14.4 

7.2 

.2 

61.3  1 

14.6 

7.3 

.0 

61.3  i 

14.6 

7.3 

-.2  | 

!  61.4  1 

14.6 

7.3 

-.4  ! 

61.4  1 

14.6 

7.3 

-.8 

60.9 

14.3 

7.1 

The  actual  X,  values  are  quite  insensitive  to  z  except  at  c 


RMS(EST)/SD[X(1)] 
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Interestingly  enough,  the  largest  effect  of  correlation  is  on  the 
light-tailed  distributions,  most  pronouncedly  for  large  positive  correlation. 

On  examining  the  statistics,  the  reason  for  the  loss  of  efficiency  is  not 
that  the  bias  has  increased.  In  fact,  the  exponential  tail  estimates  have 
very  little  bias  at  p  =  +  .8  for  the  light-tailed  distributions.  The 
problem  is  that  their  variances  increase  considerably. 

At  any  rate,  the  exponential  tail  estimates  hold  up  fairly  well, 
always  have  an  RMS  error  less  than  that  of  y  and  are  a  considerable  improve¬ 
ment  uniformly  for  the  heavy  tailed  distributions.  Of  course,  exponentially 
decreasing  auto-correlation  in  our  example  rules  out  long  term  dependence, 
and  one  could  manufacture  examples  of  stationary  sequences  with  long  term 
dependence  where  the  exponential  tail  estimate  would  give  very  poor  results. 
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Appendix  I 

MEAN  AND  VARIANCE  CALCULATIONS 

This  appendix  contains  the  calculations  leading  to  the  mean  and  vari 
ance,  under  the  quadratic  model,  of  the  statistic 

m-1' 

9  =  2  kuik[X(k)  '  X(k+1)] 

1 

Put  Wq  =  0  andyk=  kwk  -  (k  -  l)^k_1  to  get 

m-1^ 

9  =  S  7kL^E(k)^ 

1 

m-1  m-1 

=  a  SykE(k)  +  \  2  YkE(k) 

1  1 

We  use  repeatedly  the  fact  that  the  E^  have  the  representation 

m-1  7 

E(k)  ■  53  -f 

J=k 


where  Zm  ^  are  independent  exponential  variables  with  mean  one. 

Therefore 


E  ^E(k)  ^  yk 


and 
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m-1 


E(E^)  =  e£ 


3  ,£=k 


ZkZ* 

~3T“ 


^4  x 

j* 


J  ,i=k 


Hence 


m-1 


2  .  (2) 

-  yk  +  uk 


m-1 


E9  = 


Yk 


For  any  sequence  8^,  k  =  1,  ....  m-1 


m- 1 


m-1 


2  Vk  =  2  K(3k 


using  the  convention  S„  =  0.  Hence 
3  m 


m-1  m-1 

^  Ykyk  =  2 


and 
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rc-1  m-1 

Syk[Mk +  =  -  4+] + 

i 


k+l  k2 


m-l 


■ 2  £  Vi, 


so 


rc-1  m- 1 

E0  =  a  22  +  b  ^  u 

1  1 


"k^k 


To  get  the  variance  of  9,  write 


HH  m-l  /m-l  ,  7 

b  / t — *  lil 

k^k  T  ? 

1  k=l  \j ,2=k 


•  •  •  i:  ■*  *  $  z>k(r 


j  r 


and  denoting 


min(j ,£) 


j£ 


gives 


m-l  m-l 


9*a  S“kzk  +  t,Z^  VT 

k*l  j,«=l 


Let  V  2k  - 


1 .  Then 
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94 


h 


l 


m- 1 

E 


j=i 


we  have 


Var(9)  =  ^  (au)k  +  2bhk)2  +  4b  S  (aaJk  +  2bhk)hkk  +  4b 

a>("k  *  hkk)l2  +  41)2  z  *4  • 

£<k 

Note  that 

hk  +  hkk  =  hk  +  ZVk 


To  compute: 


4 


tn-1 

E 

z=i 


(2)  2 

uc 


where 


y 


(2) 

l 


t 


So  finally  we  have  the  equation  for  the  variance  given  in  the  text 

m-1 

Var(0)  =  2^(au)k  +  byk  +  bcokuk)2  +  ^  uk2)'^2 

k=l  k=i 


$ 
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Appendix  II 

DERIVATION  OF  MINIMUM  VARIANCE  UNBIASED  ESTIMATORS 
AND  MINIMUM  SQUARED  ERROR  ESTIMATOR? 


In  this  appendix  the  minimum  variance  unbiased  estimators  and  the 
minimum  squared  error  estimators  are  derived.  Writing  the  variance  as 


(II. 1) 


The  conditions  for  unbiasedness,  if 


0  =  La  +  Mb 


are 


E“ksL>  Swkyk  =  M  *  ( 

The  minimum  variance  unbiased  estimators,  minimized  at  h,  are  gotten  by 
selecting  the  to  minimize  Eq.  (II.l)  for  a  given  value  of  h  under  the 
constraints  of  Eq. (II.2). 

For  any  set  of  coefficients  the  bias  at  h  is 

8  =  a(L  +  hM  -  -  h 

To  get  the  minimum  squared  error  estimate  at  h,  minimize 
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An  exact  solution  seems  formidable.  We  get  an  "almost"  solution  by 
noting  that  the  second  term  in  Eq.(II.l)  is  usually  small  in  size  compared  with 
the  first  (except,  perhaps,  for  negative  h).  Hence,  we  throw  it  away  and 
work  on  minimizing 


m-1 

E  c“k  +  h(Yk  +  ‘Vk*!2 
1 

2  2 

with  the  constraints  of  Eq.  (II.2)  or  with  B  /a  added. 
Let 

Zk=*k  +  h(Yk  +  “kV  ; 

then 


Verifying  that 


m-1 


m-1 


E  yk =  E 1 E  =  E 


E^k +  2h  E  wkwk 


gives 
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Define  the  sequence  ak>  k  *  1,  ....  m-1  to  satisfy  the  Identity 

m-1  m-1 

S“kZk  "S'Vk 

1  1 


Then 


a  ^  Zk  -  2h  ^  Zkak 


The  minimum  variance  unbiased  problem  reduces  to  Eq. (II .4): 


Minimize  ^  Zk2  subject  to 


Z.  =  L  +  2hM 


Z“kzk  ■ M 


(II. 4) 


The  minimum  squared  error  problem  becomes  Eq . ( I I . 5): 


Minimize 


(E^i2) 


where 


TOO 


B  =^L  +  hM  Zk  +  h^akZk  j  .  (II. 5) 

In  the  first  problem,  the  solution  is  clear 

zk  =  XQ  +  xlak 

where  Xg,  X-j  are  determined  by 

(m  -  1)Xq  +  X^  ak  =  L  +  2hM 

»oZ“k*  X1  £“k  ’  M  '  (n-6) 


In  the  second  problem,  we  get  the  minimizing  equation 

Zk  =  ^  5 

so  the  solution  has  the  same  form  as  the  solution  to  the  first  problem.  We 

% 

can  solve  for  B  by  using 


£  “A -'(Zv'S  «v) 
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so 


B  =  L  +  hM  -  ^  Zk  +  h  2  akZk 

\ 

.  L  +  hM- -  hLak) 

gf  =  _ L  +  hM _ 

m  "  2h  Sak  +  h2  Sak2 


=  L  +  hM 

1  +2  (1  -  hak)2 

In  both  cases  the  solution  hinges  on  solving  forehand  then  being  able  to 
use  the  Zk  values  to  get  the  values. 

> 

l 

Write  ^  ojk  so 

1 

I  <*>k  =  Ask  =  sfe  -  .  sQ  =  0 


> 


Zl,  =  (1  +  hy.  )As. 


+  hsk 
+  IT 


» 


Then 
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and  the  identity 


Write,  assuming 


define  A+8k  =  3k 


The  identity  will 


or 


Denote  qk  =  1  + 


£Vk  =23  can  wntten  as 

Vi  ■ 


rn-1  m-1 

12  VkAsk  *  12  sk(vi  ■  “ktivi* 

l  i 

-  Sk+i  ,  so  getting  the  equation 


" 

A+  [  (1  +  huk)ak]  + 


follow  if 


hak  i 

A+CO  +  ^k^k^  +  ~T~  =  k" >  ^  =  ^  > 


m-1 


KA-C(1  +  hyk)ak]  +  hak  =  1 »  k  =  1 »  • • • »  m_1 
uk  and  look  at  the  homogeneous  equation  (8  =  0) 

KA+(qkSk)  +  h8k  =  0,  k  =  1 ,  m-1 


(II. 


The  solution  B, 


or 


Since  uj^  +  As^ 
Bk  sequence  is 

xk*  Yk 


Thus 


to  this  equation  has  two  uses:  First,  write 


hB, 

“IT 


] 


hl 

F 


1 


(I 


,  Eq.  (II. 8)  gives  aik  in  terms  of  Z^.  The  other  use  of  the 
this:  Let  ak  =  ak3k ,  then  note  that  for  any  two  sequences 


A+XkYk  =  XkYk  ‘  Xk+lYk+l 

-  Wk +  WA  • 


kA+(qk°kek)  +  hakek 

-  kakA-qk6k  +  kVl6k+lA-ak  +  hak6k 
•  kqk+i3k+iVk  • 


To  solve  Eq. (II. 7),  then,  we  need  to  solve 
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k  Vl  Bk+1  A+ 


=  1 


or 


A+  ak 


1 

kqk+l  Sk+1 


To  solve  for  3k>  write 


A+  (qk3k)  ■+  ^  =  0 


or 


or  for  q  k  >  1 , 


^qk  +  k^k  =  qk+l  8k+l 


k-1 

n  (q, 


+  £) 

J  J 


Define  nk,  k  =  0, 


m-1  by  Hg  =  1 


\  ■  ?  (1  +  35J> 


Since  6-|  is  arbitrary,  use  the  solution 
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r 


r 


» 


t 


t 


» 


Then 

8»<"t  +  I>  ■  "*-l  <T^>  =  h  • 

so  Eq.  (I I. 8)  becomes 


Going  back  to  the  ak  sequence,  we  have 


Vk  = 


1 


kg. 


k+1  4k+l 


kn, 


leading  to  the  solution 


and  finally,  to. 


thus  completing  the  work. 

The  03 k  given  by  the  above  are  not  easy  to  compute  by  hand;  a 
short  computer  program  was  written  to  calculate  them.  The  question  then 
came  up  of  whether  the  simple  approximate  solution  given  in  Subsection  3.2. 
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1 


was  close  enough  to  the  complicated  solution  given  above  so  that 
be  used  with  almost  equal  effectiveness. 

The  first  criterion  we  computed  was  average  percentage  of 
i.e.,  i f  o> ,  is  the  above  solution  at  h  andui,'  is  the  approximate 

In  < 


if  could 

di fference ; 
solution , 


Avg.  %  Diff.  =  100  x 

This  is  tabulated  below. 


Average  Percentage  of  Difference 


h  = 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

.0 

-0.1 

%  Diff. 

14.6 

10.9 

10.5 

9.9 

8.9 

7.4 

4.9 

.  00 

-14.6 

The  approximation  and  the  solution  given  above  diverge  for  negati  /e 
h.  This  may  be  due  to  the  fact  that  the  second  term  of  Var(r),  which  we 
discarded,  becomes  significant  for  negative  h. 

To  gauge  the  effects  of  this  difference  we  ran  both  solutions  at 
h  =  0.5  and  h  =  0.25  as  estimators  of  X^.  The  comparison  is  given  below. 


Comparison  of  Ratios 


b  = 

T~ 

.75 

1 .0 

1.25 

1 .50 

1.75 

2.0 

Exact  .5 

.46 

1.12 

2.47 

3.69 

4.73 

5.60 

6.35 

Weibull 

Approx  .5 

.44 

1 .09 

2.33 

3.46 

4.42 

5.21 

1 

5.90 

Weibul 1 

Exact  .25 

.64 

.56 

1 .39 

2.23 

2.96 

3.56  i 

4.10 

Weibul  1 

Approx  .25 

.61 

.56 

1.34 

2.12 

2.80 

3.36  j 

3.85 

Weibul! 

Exact  .5 

.54 

.37 

.42 

.88 

1 .37 

1 .81  ! 

2.21 

Lognormal 

Approx  .5 

.52 

.36 

.41 

.84 

i 

1 .30 

j 

1 .71 

2.07 

Lognormal 

Exact  .25 

.59 

.52 

.35 

.43 

.70 

.98 

1 .24 

Lognormal 

Approx  .25 

.56 

.50 

i 

.34 

.43 

.68  | 

.94  | 

1.19 

! 

Lognormal 

These  are  close,  and  the  biases  are  equally  close.  Thus,  for  long¬ 
tailed  distributions,  which  is  the  only  solution  in  which  the  linear 
model  produces  effective  estimators,  the  approximate  solution  probably 
can  be  substituted  for  the  exact  solution  without  a  y  deterioration  in 
performance. 


t 
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Comparison  of  Ratios 


b  = 

.5 

.75 

1.0 

1 .25 

1  .50 

1.75 

2.0 

Exact  .5 

.46 

1.12 

2.47 

3.69 

4.73 

5.60 

6.35 

Weibul  1 

Approx  .5 

.44 

1 .09 

2.33 

3.46 

4.42 

5.21 

5.90 

Weibul 1 

Exact  .25 

.64 

.56 

1 .39  | 

2.23 

2.96 

3.56 

4.10 

Weibul  1 

Approx  .25 

.61 

.56 

1.34 

| 

2.12 

2.80 

3.36 

3.85 

Weibul! 

Exact  ,5 

.54 

.37 

.42 

.88 

1 .37 

1 .81  i 

2.21 

Lognormal 

Approx  .5 

.52 

.36 

.41 

.84 

1 .30 

1.71 

i 

2.07 

Lognormal 

Exact  .25 

.59 

.52 

.35 

.43 

.70 

.98 

1  .24 

Lognormal 

Approx  .25 

.56 

i _ 

.50 

.34 

i _ 

.43 

.68 

lJ!_ 

1.19 

Lognormal 

These  are  close,  and  the  biases  are  equally  close.  Thus,  for  lonq- 
tailed  distributions,  which  is  the  only  solution  in  which  the  linear 
model  produces  effective  estimators,  the  approximate  solution  probably 
can  be  substituted  for  the  exact  solution  without  any  deterioration  in 
performance. 


